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ABSTRACT 


Algebraic  multigrid  (AMG)  is  a  numerical  method  used  to  solve  particular 
algebraic  systems,  and  interest  in  it  has  risen  because  of  its  multigrid-like  efficiency. 
Variations  in  methodology  during  the  interpolation  phase  result  in  differing  conver¬ 
gence  rates.  We  have  found  that  regular  interpolation  weight  definitions  are  inade¬ 
quate  for  solving  certain  discretized  systems  so  an  iterative  approach  to  determine 
the  weights  will  prove  useful.  This  iterative  weight  definition  must  balance  the  re¬ 
quirement  of  keeping  the  interpolatory  set  of  points  “small”  in  order  to  reduce  solver 
complexity  while  maintaining  accurate  interpolation  to  correctly  represent  the  coaxse- 
grid  function  on  the  fine  grid.  Furthermore,  the  weight  definition  process  must  be 
efficient  enough  to  reduce  setup  phase  costs. 

We  present  systems  involving  matrices  where  this  iterative  method  signifi¬ 
cantly  outperforms  regular  AMG  weight  definitions.  Experimental  results  show  that 
the  iterative  weight  definition  does  not  improve  the  convergence  rate  over  standard 
AMG  when  applied  to  M-matrices;  however,  the  improvement  becomes  significant 
when  solving  certain  types  of  complicated,  non-standard  algebraic  equations  gener¬ 
ated  by  irregular  operators. 
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I. 


INTRODUCTION 


Precursors  to  multigrid  (MG)  methods  were  developed  in  the  mid  1960’s  by 
Fedorenko  (1964),  and  Bachvalov  (1966),  and  remained  fairly  unknown  until  Achi 
Brandt  [Ref.  1]  fully  developed  multigrid  in  his  1973  article.  In  that  paper,  he 
demonstrated  for  the  first  time  the  utility  of  true  multigrid  methods  as  a  fast  numeri¬ 
cal  solvers  for  boundary  value  problems.  An  explosion  of  papers  ensued  and  the  1980’s 
produced  a  plethora  of  articles  as  well  as  efficient  and  reliable  computer  algorithms  to 
further  demonstrate  the  practicality  and  usefulness  of  MG.  Confidence  in  multigrid 
grew  when  the  true  merits  of  this  solver  came  to  bear  fruit  in  faster  convergence  rates 
over  conventional  methods.  Proofs  of  these  faster  convergence  rates  coupled  with  a 
sufficient  number  of  satisfactory  numerical  results  gave  way  to  a  general  acceptance 
of  the  method  even  to  the  most  sceptic  [Ref.  2]. 

Originally  developed  to  solve  simple  boundary  value  problems,  multigrid  meth¬ 
ods  gained  wide  recognition  for  their  speed  and  efficiency  in  solving  general  linear 
partial  differential  equations  (PDFs)  of  the  elliptic  form,  e.g.,  (in  two  dimensions) 

-V\{x,y)  =  f{x,y), 

or  explicitly  as 

Uxx  '^yy  ~ 

If  this  equation  is  now  discretized,  for  example  by  finite  differences,  on  some  rectan¬ 
gular  domain,  C  7?.^  as  in  Figure  1,  and  if  we  denote  the  boundary  of  the  domain 
as  dCl,  then  multigrid  will  solve  problems  involving  Laplace’s  equation 


— =  0  in  0  1 


u  =  0  on  do, 


(1.1) 


1 


Figure  1.  Rectangular  domain  for  discretization  of  a  2-D  problem. 


Poisson’s  equation 

=  s  in  n' 

u  =  0  on  oCl 

Helmholtz’s  equation 

— +  K^u  =  5  in 
u  =  0  on  do, 

or  even  the  anisotropic  Poisson  equation 


d^u  d^u  . 
u  =  0  on  50 


(1.2) 


(1.3) 


(1.4) 


in  any  number  of  dimensions  almost  effortlessly  (comparatively  speaking,  of  course). 

These  types  of  equations  frequently  arise  in  physical  applications  such  as 
steady  state  temperature  problems,  fluid  flow,  orbital  mechanics  (gravitational  fields), 
steady  state  electrical  field  problems,  and  many  others.  When  a  matrix  equation  arises 
from  the  discretization  of  (I.l)  through  (1.4),  multigrid  methods  are  quite  successful. 

MG  methods  are  fast  iterative  solvers  using  a  hierarchy  of  levels,  or  grids,  and 
may  be  used  with  many  types  of  discretization  techniques  such  as  finite  difference 
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or  finite  element  methods.  In  solving  problems  of  the  elliptic  nature  using  these 
discretization  techniques,  multigrid  has  proven  to  be  the  fastest  numerical  solution 
technique  in  the  field.  Unlike  other  numerical  solvers,  multigrid  is  general  enough  so 
that  it  can  effectively  use  fairly  arbitrary  regions  and  boundary  conditions  and  more 
importantly,  does  not  depend  on  the  separability  of  the  differential  equations. 

In  fact,  from  presentations  at  the  Copper  Mountain  Conference  on  Multigrid 
Methods  and  from  various  papers  which  can  be  found  on  MG  Net  (“http://casper.cs. 
yale.edu/mgnet/www/mgnet-papers.html#  Y”),  we  know  that  multigrid  can  also  be 
directly  applied  to  more  complicated,  non-symmetric  and  nonlinear  systems  of  equa¬ 
tions,  like  the  Lame-System  of  elasticity  or  the  Stokes  (or  Navier-Stokes)  equations. 
Multigrid  has  been  applied  successfully  to  electrostatic  and  magnetostatic  problems 
[Ref.  3],  to  statistical  physics  problems,  to  integral  equations  and  to  image  recon¬ 
struction  algorithms  [Ref.  4,  5].  It  can  also  be  applied  to  problems  in  control  theory, 
partical  physics  and  permeable  magnetic  materials  [Ref.  2]. 

The  main  concept  of  multigrid  methods  is  “to  complement  the  local  exchange 
of  information  in  point- wise  iterative  methods  (on  each  level)  by  a  global  one  utilizing 
several  related  systems,  called  course  levels,  with  a  smaller  number  of  variables  [Ref. 
6].”  The  key  to  multigrid’s  performance  then  is  to  apply  a  relaxation  technique 
as  many  times  as  needed  to  dampen  the  oscillatory  modes  of  the  error  (since  we 
know  relaxation  works  best  on  highly  oscillatory  modes)  and  then  apply  a  coarse-grid 
correction  scheme  to  eliminate  the  smooth  components  of  the  error.  This  combination 
of  relaxation  and  coarse-grid  correction  is  the  essence  of  why  MG  works  so  well. 

Using  multigrid  methods,  we  can  solve  a  large,  difficult  problem  by  reduc¬ 
ing  it  iteratively  into  successive  smaller,  easier  ones.  At  the  coarsest  level  then  the 
system  can  be  solved  directly  using  a  known,  efficient  direct  method  (such  as  the 
LU- decomposition,  Gaussian  elimination,  etc.).  For  systems  such  as  (I.l)  through 
(1.4),  very  efficient  methods  are  available.  From  A  Multigrid  Tutorial  [Ref.  7],  we 
know  that  when  applied  to  am  N  x  N  grid,  multigrid  methods  are  nearly  optimal 
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since  they  only  require  0{N^  log  N)  arithmetic  operations  and  hence  approach  the 
minimum  operation  count  of  0{N'^)  operations.  In  fact,  it  is  shown  in  [Ref.  7]  that 
the  full  multigrid  V-cycle  (FMV)  method  which  we  will  discuss  in  the  next  chapter, 
requires  only  0{N^)  to  the  level  of  discretization  on  the  “model”  problem  (Laplace’s 
equation  (LI))  of  dimension  d.  This  operation  count  is  optimal. 


II.  THE  FUNDAMENTALS  OF 
MULTIGRID 


There  are  two  fundamental  components  of  multigrid;  one  is  the  idea  of  coarse- 
grid  correction  (CG)  and  the  other  is  that  of  nested  iteration  (NI).  Presented  here  is  a 
basic  outline  of  both  methods.  To  begin,  we  note  that  the  principle  concept  of  CG  is 
that  it  is  a  correction  strategy  that  transfers  the  components  of  the  problem  between 
the  fine  and  coarse  grids.  Complementing  coarse-grid  correction  is  nested  iteration 
which  is  based  on  the  following  concept:  use  the  information  on  the  coarse  grids  to 
provide  an  informed  guess  for  the  initial  guess  on  the  finer  level. 


A.  THE  PROBLEM  STATEMENT 

In  order  to  understand  the  basics  of  multigrid  and  to  motivate  its  usage,  let 
us  apply  the  method  to  a  simple  one-dimensional  problem.  We  begin  by  discretizing 
a  problem  of  the  form 


-  +  (^u{x)  =  f{x),  0  <x  <1,  a>0  (III) 

with  Dirichlet  boundary  conditions  {u{x)  =  0  on  dQ)  by  partitioning  the  domain,  as 
in  Figure  2,  into  N  +  1  points: 


Xj  =  jh,  where  j  =  0,...,N. 


We  make  h  —  1/N  of  constant  width  throughout  the  interval.  This  creates  a  grid 


x=0 

G - e - © - e- . —  o - ©- 

Xo  Xi  X2  Xj 


X=1 

-e - o 

%-i  % 


Figure  2.  Discretized  domain  of  problem  ILL 


of  N  —  1  interior  points. 


At  each  of  these  interior  points,  the  equation  in  (II.l)  is 
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replaced  by  a  second  order  discretization,  such  as  a  finite  difference  approximation 

(11.2) 

Vo  =  =  0 


where  we  define  fj  =  f{xj).  The  approximation  (II.2)  of  (II.l)  leads  to  a  system 
of  linear  equations.  In  compact  form,  it  is  written  as  Av  =  f,  or  in  explicit  matrix 
notation,  it  is  written 


2 -f  <7^^  —1 

Vl 

h 

1 

-1  2  + ah?  -1 

V2 

f2 

*  •  ’  •  ’  . 

; 

• 

h? 

-1  2  A  ah?  -1 

-1  2  -h  ah? 

ViV-1 

/n-i 

where  the  matrix  A  is  tridiagonal,  symmetric,  positive  definite  and  has  dimension 
N-lxN-1.  u  is  a  coltunn  vector  containing  the  approximate  discretized  solution 
of  u{x)  and  /,•  =  f(xi)  is  the  discretized  right-hand  side  (RHS). 


B.  NOTATION 

Before  we  continue  with  the  outline  it  is  instructive  to  define  a  few  terms.  Let 
denote  the  original  grid  or  the  domain  on  which  the  original  problem  is  defined. 
We  call  the  finest  grid.  The  symbol,  is  defined  to  mean  the  next  coarser  grid 
(usually  of  size  iV/2  —  1  but  certainly  not  restricted  to  that).  Similarly,  u*  will  be 
an  approximation  to  the  exact  solution,  on  Qfi.  This  means  that  represents 
the  approximate  solution  to  on  the  next  coarser  level  This  notation  will  be 
consistent  for  all  such  grids,  The  algebraic  error,  on  (  A:  =  2”,  n  = 
0, 1, 2, 3, . . .)  is  given  then  by  where  A;  =  1  represents  the  finest  level. 

When  we  translate  the  vectors  only  between  two  consecutive  grids,  we  will  use 
the  notation  to  represent  the  finer  of  the  two  grids  and  Qfi  to  denote  the  coarser 
one.  This  maLes  it  easier  to  follow  the  grid  transfer  sequence  vice  getting  caught  up 
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in  exponential  notation  of  the  varying  grids.  Now  that  we  have  dispensed  with  some 
of  the  notational  formalities,  let  us  move  on  to  more  of  the  important  issues. 

C.  ITERATIVE  METHODS 

We  introduce  iterative  methods  here  to  illustrate  the  need  for  coarse-grid  cor¬ 
rection.  First  of  all,  iterative  methods  are  needed  and  in  fact,  are  required  since  the 
use  of  direct  methods  can  be  impractical  if  the  matrix  A  is  large  and  spaxse.  The 
sought-after  factors  of  A  using  these  direct  methods  can,  and  tend  to,  be  very  dense. 
Even  when  the  matrix  is  banded  (but  still  sparse),  algorithms  that  are  ideal  for  fac¬ 
toring  such  problems  may  be  difficult  to  implement  [Ref.  8].  Iterative  methods  by 
contrast  “generate  a  sequence  of  approximate  solutions  and  essentially  involve  the 
matrix  A  only  in  the  context  of  matrix- vector  multiplication  [Ref.  8].”  Iterative 
methods  are  attractive  and  easy  to  use  because  of  their  simplicity  in  implementation 
but  may  be  prohibitively  slow  when  the  error  is  smooth. 

What  does  it  mean  when  we  say,  “when  the  error  is  smooth?”  To  answer  this 
and  to  understand  the  need  for  coarse-grid  correction,  we  must  first  show  the  power 
of  iterative  methods  and  then  focus  on  their  limitations.  To  avoid  any  confusion,  we 
note  here  that  the  phrase  “iterative  method”  is  isomorphic  to  “relaxation  scheme” 
or  just  simply,  “relaxation.”  We  demonstrate  relaxation  with  the  weighted  Jacobi 
iteration  since  it  is  simple  in  nature  and  offers  the  same  benefits  as  other  iterative 
methods;  in  addition,  it  shows  the  common  defects  of  relaxation  in  general.  Other 
iterative  methods  include  (regular)  Jacobi,  Gauss-Seidel,  successive  over-relaxation 
(SOR)  and  Chebyshev  semi-iterative;  of  course,  there  are  many  others. 
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1.  The  Weighted  Jacobi  Iteration  on  the  “Model” 
Problem 

To  begin,  let  us  solve  the  following  one-dimensional  “model”  problem  which 
is  problem  (II-2)  with  a  and  /  set  to  zero.  The  problem  then  becomes 


-Uj-i  +  2uj  —  Uj+i  =0,  1  <  j  <  N  -1 

Uq  =  U]^  =  0. 


(11.3) 


We  solve  (II.3)  by  using  the  weighted  Jacobi  iteration.  This- iteration  is  computed 
using 

where 

=  I  ,  l<j<N-l 

and  a;  €  is  a  weighting  factor  chosen  such  that  0  <  a;  <  1.  In  our  use  of  the 
weighted  Jacobi  iteration,  we  set  oj  =  2/3.  It  turns  out  to  be  the  optimal  weighting 
factor. 

For  the  moment,  let  us  assume  that  we  have  found  an  approximation,  v  ^ 
0,  to  the  discretized  solution,  u  =  0,  on  our  model  problem.  Since  the  system  is 
homogeneous  and  linear,  we  know  the  error,  e,  exactly,  namely  —v  since  e  =  u  —  v  = 
—V.  In  reality,  the  error  will  consist  of  many  different  frequencies  but  for  simplicity, 
let  us  assume  that  e  consists  of  only  three  frequency  modes:  one  high,  one  low  and  a 
third  mode  in  the  middle.  Our  assumption  using  the  three  modes  here  is  merely  to 
amplify  the  benefits  of  relaxation  while  simultaneously  exploiting  the  inherent  defects. 

2.  Frequency  Modes  of  the  Error 

Consider  an  initial  approximation  to  the  solution  consisting  of  the  Fourier 

modes 


where  0  <  J  ^  TV  and  1  <  A:  <  TV  —  1.  The  number  k  represents  the  wave  number 
of  the  Fourier  mode.  The  fcth  mode  consists  oi  kj2  full  sine  waves  on  the  interval. 
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As  we  will  see  shortly,  the  wave  number  k  (frequency  mode)  will  play  a  major  role 
in  the  overall  effect  of  relaxation.  Figure  3  illustrates  the  modes  Vj  =  sin  with 
wave-numbers  k  —  6,  12  on  a  grid  with  N  =  64.  Wave  number  k  =  1  represents 

the  lowest  frequency  while  A:  =  63  is  the  highest  one.  We  now  iterate  (or  “relax”)  on 


Figure  3.  Fourier  modes  Vj  =  sin  with  wave-numbers  k  =  1,  6,  12  ore  a  grid 
with  N  =  64.  The  kth  mode  consists  of  kj2  full  sine  waves  on  the  interval.  vi  is 
represented  by  the  solid  line,  vq,  the  dashed  and  Vi2,  the  dot-dashed. 

problem  (II.3)  by  applying  the  weighted  Jacobi  iteration  with  to  =  2/3  on  a  grid  of 
size  N  =  64. 

But  wait!  Maybe  we  axe  being  a  bit  too  hasty.  Let’s  step  back  from  solving  the 
problem  for  a  moment.  Before  we  see  what  happens  with  our  mixed-frequency  initial 
guess,  it  is  instructive  to  first  watch  what  happens  to  the  individual  components  in 
our  initial  guess.  Figure  4  shows  the  application  of  the  weighted  Jacobi  iteration 
metho  on  problem  (II.3)  with  Vi,  Vz,  Vq,  and  Vi2  as  separate  initial  guesses. 
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Iterations 


Figure  4.  The  weighted  Jacobi  iteration  with  a;  =  2/3  applied  to  the  1-D  ^^modeF 
problem  with  iV  =  64  and  with  initial  guesses  Ui,  vz,  vq,  Vi2.  ||el|oo  is  plotted  against 
the  iteration  number  for  50  iterations. 


Notice  that  the  error  quickly  decreases  to  zero  within  a  few  iterations  for  the 
high-  frequency  mode  (^12)  and  although  not  quite  as  rapidly,  the  error  also  decreases 
for  the  middle-frequency  mode  (ue)  but  it  doesn’t  quite  go  to  zero.  In  extreme  con¬ 
trast,  the  low-frequency  modes  (ui  and  U3)  appear  to  be  relatively  untouched  by  the 
relaxation  scheme.  In  fact,  a  prohibitively  large  number  of  iterations  will  still  not 
effectively  reduce  the  error.  It  is  the  wonderful  property  of  iterative  methods  which 
allow  for  the  quick  elimination  of  high-frequency  modes.  On  the  other  hand,  they  are 
quite  defective  in  decreasing  the  error  when  it  consists  solely  of  low-frequency  modes. 

We  now  return  to  problem  (II.3)  with  our  mixed-frequency  initial  guess.  It  is 


of  the  form 


-h  sin 


-f  sin 
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Figure  5  shows  how  the  error  of  a  typical  approximation  vector  decreases  using  re¬ 
laxation  alone.  It  is  not  uncommon  that  the  error  decreases  rapidly,  usually  within 
a  few  iterations.  This  is  due  to  the  quick  elimination  of  the  high-frequency  modes. 
But  after  this  initial  rapid  decrease,  the  error  is  reduced  much  more  slowly.  The 
main  culpritsare  the  low-frequency  modes.  “The  important  observation  is  that  the 
standard  iterations  converge  very  quickly  as  long  as  the  error  has  high-frequency  com¬ 
ponents.  However,  the  slower  elimination  of  the  low-frequency  components  degrade 
the  performance  of  the  methods  [Ref.  7].” 


Figure  5.  The  weighted  Jacobi  iteration  with  uo  =  2/3  applied  to  problem  (ITS)  with 
iV  =  64  and  with  the  initial  guess  Vj  =  |  [sin  +  sin  +  sin  (^)]  •  ||e||oo  is 
plotted  against  the  iteration  number  for  50  iterations. 

What  we  have  seen  is  that  iterative  methods  work  very  well  for  the  first  several 
iterations  but  inevitably  the  convergence  rate  slows  and  the  methods  seem  to  stall. 
This  phenomenon  is  due  to  the  fact  that  the  rapid  decrease  in  error  in  the  early 
stages  of  the  method  is  from  the  swift,  efficient  elimination  of  the  oscillatory  modes 
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(high-frequency  components).  However,  once  the  high-frequency  components  have 
been  removed  from  the  error,  the  iteration  is  less  effective  in  reducing  the  remaining 
smooth  modes  (low-frequency  components). 

This  property  of  quickly  removing  the  oscillatory  modes  but  leaving  the  smooth 
modes  of  the  error  is  called  the  smoothing  property.  The  smoothing  property  of 
iterative  methods  (relaxation)  is  a  serious  limitation.  How  do  we  overcome  this  grave 
obstacle?  It  is  in  the  application  of  coarse-grid  correction  to  the  remaining  “smooth 
error”  that  we  remedy  this  problem. 

D.  THE  METHOD  OF  COARSE-GRID  CORRECTION 

In  light  of  the  limitations  of  iterative  methods,  we  seek  a  post-relaxation 
scheme.  The  purpose  of  the  coarse-grid  correction  scheme  (CG)  is  to  eliminate  the 
low-frequency  modes  of  the  error  once  relaxation  has  become  ineffective.  The  means 
by  which  we  do  this  is  intergrid  transfers.  Intergrid  transfers  move  vectors  (and  the 
matrix)  from  the  fine  grid  to  the  coarse  grid  and  vice-versa  by  means  of  restriction 
and  interpolation  operators.  In  “moving”  to  a  coarser  grid,  we  observe  that  the 
smooth  modes  of  the  error  become  higher-frequency  modes  on  this  new  grid.  On  the 
coarse-grid  now,  the  error  has  oscillatory  components  again;  but  we  can  effectively 
remove  those  modes  using  relaxation.  This  two  step  methodology  is  the  basis  for  CG. 
To  understand  what  it  means  to  “move”  a  vector  between  grids,  we  introduce  a  few 
required  operators. 

1.  The  Restriction  Operator 

Starting  on  the  fine  grid,  we  move  the  problem  to  the  next 

coarser  grid,  by  applying  an  intergrid  transfer  to  A^,  and  /^.  This  operation 
gives  us  the  coarse-grid  equivalent  to  that  problem  but  on  In  doing 

this  though,  we  must  be  careful  to  ensure  that  the  coarse-grid  problem  “be  consistent 
with  the  differential  problem  in  the  same  way  as  the  fine-grid  problem  [Ref.  2].”  The 
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action  of  moving  to  the  coarse-grid  problem  is  performed  by  the  restriction  operator 

/f 

is  a  linear  operator  from  to  and  has  rank  iV/2— 1  (see  [Ref.  2,  7, 9]  for 
more  on  intergrid  transfers).  The  most  obvious  of  these  types  of  restriction  operators 
is  injection.  Injection  is  a  linear  operator  which  produces  the  coarse-grid  solution 
approximation 

=  /f 

where 

vf  =  v^j  for  l<i<Y-l. 

In  words,  injection  gets  the  value  of  the  coarse-grid  vector  directly  from  the  associated 
fine  grid  point.  Figure  6  shows  how  the  injection  operator  acts  on  a  discretized  sine 
wave  when  moving  the  vector  from  to 

Another  type  of  restriction  operator  is  full  weighting.  It  is  also  a  linear  operator 
from  to  with  rank  N/2  —  1.  Like  injection,  it  produces  the  coaxse-grid 

solution  approximation 

but  in  this  case  the  transfer  to  the  coarse  grid  is  a  bit  different: 

=  \  H'-i  +  +  4i+i)  for  1  <  j  <  Y  -  1. 

In  the  case  of  the  full  weighting  restriction  operator,  “the  values  of  the  coarse  grid 
vectors  are  a  weighted  average  of  values  at  neighboring  fine  grid  points  [Ref.  7].” 
Figure  7  shows  how  the  full  weighting  restriction  operator  acts  upon  a  vector  when 
moving  the  vector  from  the  fine  grid  to  the  coarse  grid. 
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Figure  6.  Restriction  hy  injection  on  a  discretized  sine  wave  from  to  .  The 
circles  (o)  on  Of  directly  become  the  coarse-grid  points  on  . 

2.  The  Interpolation  Operator 

Once  on  the  coarse  grid,  and  after  we  complete  our  calculations  there,  we 
need  to  be  able  to  return  to  We  introduce  the  tool  that  will  assist  us  in  this  task: 
the  interpolation  operator 

The  interpolation  operator  is  a  linear  operator  from  to  and  has  full  rank. 

Iff  moves  the  coarse-grid  vectors  up  one  level  by  the  calculation 

where 

4  =  t;f  for  0<i<y-l. 
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Figure  7.  Restriction  by  full  weighting  on  a  arbitrary  vector  from  to  .  The 
circles  (o)  on  are  the  coarse-grid  points  and  are  weighted  averages  of  neighboring 
fine  grid  point  values. 

The  interpolation  operator  produces  a  vector  on  the  fine  grid  that  is  a  weighted 
average  of  its  adjacent  points  on  the  coarse  grid.  Figure  8  graphically  depicts  the 
action  of  on  a  discretized  cosine  wave. 

Now  that  we  have  our  intergrid  transfer  operators  defined,  it  makes  sense  to 
talk  about  moving  the  matrix  between  levels.  On  fi'",  =  A.  But  what  is  A^ 

and  how  do  we  get  it?  As  a  definition,  we  will  call  A^  the  version  of  A’^  or,  the 
coarse-grid  operator.  We  get  A^  from  A^  by  the  calculation 

A^  =  (II.4) 

and  we  recover  A^  from  A^  by  the  calculation 
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Figure  8.  Interpolation  of  a  cosine  wave  from  to  The  stars  (^)  on  Of  are  the 
interpolated  values  from  their  adjacent  coarse-grid  points  on  , 

This  is  how  we  move  the  operator  A  on  all  levels.  The  coarse-grid  operator  (IL4)  is 
called  the  Galerkin  condition.  A  second  special  and  important  relationship  between 
the  full  weighting  restriction  operator  and  the  linear  interpolation  operator  is 

/^  =  c {jhf  where  cell;  (11.5) 

that  is,  they  axe  transposes  of  each  other  up  to  a  (real)  constant  which  turns  out  to 
be  essential.  (11.4)  along  with  (II.5)  are  called  the  variational  properties. 

We  note  here  that  there  is  a  disadvantage  in  using  a  linear  interpolation  op¬ 
erator.  In  fact,  Wesseling  explains: 

Because  of  the  arbitrariness  in  choosing  the  direction  of  the  diagonals  of  the 
interpolation  triangles,  there  may  be  a  loss  of  symmetry,  that  is,  if  the  ex¬ 
act  solution  of  a  problem  has  a  certain  symmetry,  it  may  happen  that  the 
numerical  solution  does  not  reproduce  this  symmetry  exactly,  but  only  with 
truncation  error  accuracy.  [Ref.  2] 
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However,  linear  interpolation  is  much  cheaper  (and  easier)  to  use  than,  for  example, 
bilinear  or  quadratic  interpolation.  So  in  our  work,  we  only  use  linear  interpolation. 

3.  The  Algorithm  of  CG 

To  begin  the  coarse-grid  correction  algorithm,  we  perform  a  relaxation  scheme 
(e.g.,  Gauss-Seidel  or  Jacobi)  on  on  the  original  equation  =  /*  with  an 

arbitrary  initial  guess  v’^.  This  has  the  effect  of  producing  a  much  better  guess 
since  it  removes  the  oscillatory  components  of  the  error  right  from  the  start.  We 
now  compute  the  residual  vector,  r^,  on  the  coarser  grid  r^  =  and 

then  solve  A^e^  =  r^  for  the  algebraic  error,  on  Cl^.  Now  that  the  error  is 
computed,  we  return  to  by  using  the  linear  interpolation  operator,  7^,  defined 
above.  Since  the  error  was  smooth  on  transfers  that  error  very  accurately 

back  to  We  now  correct  the  fine  grid  approximation  by  adding 

the  interpolated  error  e*  (=  I^e^)  back  into  the  initial  guess.  Since  A  e^, 

we  have  improved  our  fine  grid  approximation  to  the  solution  and  in  theory,  we  have 
found  the  solution.  To  improve  it  further,  we  perform  relaxation  on  the  original 
equation  A^u^  =  with  the  updated  guess  to  obtain  a  final  approximation  to  the 
exact  discretized  solution  u^.  The  coarse-grid  correction  scheme  (in  compact  form) 
is  printed  below  and  is  reproduced  from  Briggs’,  A  Multigrid  Tutorial  [Ref.  7]. 

^CG{v^,f) 

Relax  i/i  times  on  =  /*  on  with  initial  guess  v^. 

Compute  r^  =  —  A^v^). 

Solve  A^e^  =  r^  on  0^. 

Correct  the  fine  grid  approximation:  A  • 

Relax  1/2  times  on  A^u^  =  on  fi*  with  (updated)  initial  guess  v^. 

Here,  and  V2  are  small  positive  integers,  usually  between  one  and  three. 
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E.  THE  MULTIGRID  V-CYCLE  SCHEME 

The  multigrid  V-cycle  scheme  is  a  recursive  algorithm  that  iteratively  imbeds 
the  two  level  CG  algorithm  within  itself.  The  physical  property  of  this  scheme  (the 
pattern  it  creates)  is  how  it  gets  its  name.  Figure  9  shows  a  graphic  illustration  of 
the  V-cycle  pattern. 


Figure  9.  The  multigrid  V-cycle. 


We  start  one  CG  algorithm  on  the  finest  level  and  a  new  one  on  each  successive 
level  all  the  way  down  to  the  coarsest  level  where  we  the  problem  is  solved  directly. 
After  we  solve  the  problem  on  the  coarsest  level,  we  now  complete  each  CG  algorithm 
at  each  successive  level  all  the  way  back  to  the  finest  level.  To  clarify  what  the  V-cycle 
is  really  doing,  we  present  the  algorithm  in  a  more  structured  format.  In  the  code, 
it  is  assumed  again  that  there  are  Q  >  1  grids  with  the  coarsest  grid  spacing  given 
by  Lh,  where  L  =  2^~^.  This  procedure  is  also  a  reproduction  from  Briggs’  tutorial 
[Ref.  7]. 
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v^^V{v\f) 

Relax  on  ui  times  with  initial  guess 

Compute 

Relcix  on  ui  times  with  initial  guess  =  0. 

Compute 

Relax  on  =  p^  i/i  times  with  initial  guess  =  0. 

Compute  p^  =  Iiy^. 

Solve  A^^u^^  =  p^. 

Correct  •« — 

Relax  on  A*^u^^  =  p^  i/2  times  with  initial  guess 

Correct  ^ 

Relax  on  A?^u^^  =  p^  1/2  times  with  initial  guess 

Correct 

Relax  on  A^P  =  p  1/2  times  with  initial  guess  v^. 

For  all  its  simplicity,  this  algorithm  is  the  basic  structure  of  the  multigrid 
methodology.  It  takes  individual  techniques  and  well  known  concepts  (relaxation  and 
intergrid  transfers)  inherent  with  individual  defects  and  integrates  them  in  a  cohesive 
way,  capitalizing  on  the  benefits  of  each  to  produce  an  algorithm  that  is  extremely 
efficient  and  quite  powerful. 

Given  the  recursive  nature  of  the  V-cycle  algorithm  above,  we  now  present  it 
in  a  more  compact  form. 
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^  Viv\f) 

1.  Relax  i/i  times  on  with  a  given  initial  guess  v^. 

2.  If  =  the  coarsest  grid,  then  goto  4. 

Else  ^  Il\f^  - 

0; 

■y2fe  ^  y2h^y2h^  P^). 

3.  Correct  <— 

4.  Relax  V2  times  on  A^u^  =  p  with  initial  guess  v^. 

The  V-cycle  is  just  one  of  many  multigrid  cycling  schemes.  To  see  its  effec¬ 
tiveness,  we  show  in  Figure  10  how  the  V-cycle  algorithm  tailors  itself  to  the  problem 
vice  the  size  of  the  matrix  involved.  The  number  of  V-cycles  required  to  solve  the 
problem  to  truncation  error  barely  grows  even  though  we  considerably  increase  the 
matrix  size. 


Figure  10.  Size  of  the  problem  verses  number  of  V-cycles  performed. 
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F.  THE  STRATEGY  OF  NESTED  ITERATION 

When  no  information  about  the  solution  to  the  problem  is  known  before-hand, 
it  is  computationally  wasteful  to  begin  solving  the  problem  on  the  finest  grid  with  just 
any  initial  random  guess,  v^.  If  a  poor  choice  is  made,  the  speed  of  convergence  may 
be  strongly  affected.  The  algorithm  may  altogether  fail  to  converge  if  constrained 
by  the  number  of  iterations  especially  when  the  problem  is  nonlinear.  Computing 
the  residuals  on  the  “smaller”  coarser  grids  is  so  much  cheaper  and  it  makes  more 
sense  than  to  remain  only  on  O'*.  Therefore,  it  is  better  to  use  the  information  of 
the  coarse  grids  to  provide  an  informed  guess  u*  on  the  finest  grid.  In  addition  to 
providing  a  better  initial  guess  on  nested  iteration  gives  us  a  better  choice  for 
at  each  level  k  [Ref.  2].  This  idea  of  using  information  inherent  in  the  problem  from 
the  coarser  grids  to  develop  better  initial  guesses  is  the  basis  of  nested  iteration.  In 
compact  form,  the  NI  algorithm  is 

For  A;  =  2, 4, 8, 16, . . . ,  X 
Compute 

end 

Relax  1/1  times  on  Av  =  f  on  the  coarsest  level  and  solve 

J^Lh^Lh  ^  jLh  fQj. 

Interpolate  to  obtain  an  initial  guess  for  the  problem. 

For  fc  =  X  —  1, . . . ,  8, 4, 2, 1 
Compute 

Relax  1/2  times  on  Av  =  /  on  to  obtain  an  initial  guess  for  the 
problem. 

end. 

In  this  algorithm,  i/i  and  1/2  are  again  small  positive  integers  and  there  are  Q  >  1 
grids  with  the  coarsest  grid  spacing  given  by  Lh,  where  X  =  2^“^. 
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G.  THE  FULL  MULTIGRID  V-CYCLE 

In  this  section,  we  present  a  powerful  algorithm  that  combines  nested  iteration 
with  the  V-cycle.  We  also  talk  about  how  it  works  and  why  it  is  necessary.  The 
algorithm  to  which  we  are  referring  is  called  the  full  multigrid  V-cycle  (FMV)  and 
it  is  the  method  which  is  computationally  of  optimal  order.  We  already  saw  the 
algorithms  for  NI  and  the  V-cycle  and  reasoned  for  their  need.  Now  we  join  them 
into  a  single,  efficient  method. 

In  the  FMV  algorithm,  the  first  approximation  to  the  solution  is  obtained 
by  interpolation  of  the  solution  from  the  next  coarser  grid,  which  has  already  been 
computed  by  another  FMV  algorithm.  This  is  where  NI  plays  its  role.  Once  we 
achieve  a  first  approximation,  we  then  apply  the  multigrid  V-cycle.  “This  combination 
should  suffice  to  give  a  fine-grid  solution  to  the  level  of  truncation  error  [Ref.  10].” 

The  FMV  algorithm  in  compact  form  is  (reproduced  from  [Ref.  7]): 

FMV^{v^,  f) 

1.  If  =  the  coarsest  grid,  then  goto  step  3. 
else 

f2h  ^  J2h(^jh  _ 

0; 

2.  Correct  -{■  I^hP^- 

3.  u*  4—  V^{v^,  p)  vq  times. 

Figure  11  shows  how  the  FMV  scheme  operates  over  five  grids.  Notice  that 
each  V-cycle  is  proceeded  by  a  smaller  V-cycle  which  is  designed  to  provide  the  best 
initial  guess  for  the  following  cycle.  Shortly,  we  will  see  that  the  extra  work  required  in 
the  preliminary  V-cycle  is  not  only  of  minimal  cost,  but  it  will  generally  pay  for  itself 
over  the  course  in  solving  the  entire  problem  [Ref.  7].  We  also  note  here  that  FMG 
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cycles  are  less  sensitive  than  the  multigrid  cycles  but  with  one  cycle  per  refinement, 
we  can  still  guarantee  a  solution  to  the  level  of  truncation  error  [Ref.  10]. 


whereas  (o)  is  the  post  coarse-grid  correction  work.  The  (+)  on  the  first  branch 
represent  the  restriction  of  f;  no  relaxation  is  performed  here. 

Although  the  FMV  algorithm  costs  more  per  cycle  than  the  regular  V-cycle, 
it  is  also  more  effective.  “The  key  observation  in  the  FMV  argument  is  that  before 
the  problem  is  even  touched,  the  problem  has  already  been  solved  to  the  level 
of  truncation  [Ref.  7].”  The  reason  is  that  we  now  have  a  better  initial  guess  for  the 
next  finer  grid  from  our  use  of  nested  iteration.  In  fact,  we  know  that  because  of  this 
“extra”  work  during  the  preliminary  cycling  through  the  coarser  grids,  we  only  require 
0(1)  V-cycles  by  the  time  we  finally  reach  the  finest  grid.  So  the  computational  cost 
of  the  FMV  method  applied  to  a  d-dimensional  problem,  as  mentioned  before,  is 
0{N^),  which  is  optimal  vice  the  cost  of  the  V-cycle  method  alone  which  is  of  order 
OiN^^logN). 
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H.  IS  MULTIGRID  ENOUGH? 

The  heart  of  multigrid  intertwines  relaxation  methods  such  as  Jacobi  or  Gauss- 
Seidel  with  intergrid  operators  like  coarse-grid  correction.  In  consideration  of  a  gen¬ 
eral  all-purpose  elliptic  equation  solver,  the  two  techniques  enhance  each  other’s  ef¬ 
fectiveness  such  that  the  overall  algorithm  -  multigrid  -  is  more  robust  and  powerful 
than  the  sum  of  its  parts.  To  be  more  specific,  it  is  a  known  fact  that  relaxation 
effectively  smoothes  highly  oscillatory  modes  of  the  error  but  it  is  ineffective  on  the 
smooth  modes.  This  knowledge,  coupled  with  the  fact  that  any  vector  that  lies  in 
the  range  of  the  interpolation  operator  also  lies  in  the  null  space  of  the  coarse-grid 
correction  operator,  provides  us  a  process  that  ensures  that  both  the  smooth  and 
oscillatory  modes  of  the  error  are  suitably  damped. 

Multigrid  methods  have  their  limitations,  though.  Suppose  we  wish  to  solve  a 
non-elliptically  structured  problem  or  even  a  regular  elliptic  problem  but  defined  on 
an  irregular  grid.  How,  for  example,  can  we  apply  multigrid  methods  to  these  types 
of  systems.  Will  it  even  work  at  all?  Can  we  tailor  the  irregularities  in  the  meshes? 
How  do  we  apply  MG  to  problems  where  no  geometric  structiure  (or  grid)  exists? 
Some  specialized  multigrid  codes  have  been  written  for  specific,  ill-structured  prob¬ 
lems  to  satisfy  the  irregular  needs  of  the  individual  problem  but  that  defeats  the  issue 
of  providing  a  solver  that  is  more  nearly  a  “black-box”  type  of  solver.  Every  special¬ 
ization  restricts  the  scope  of  the  discretization  method  with  which  the  multigrid  code 
can  be  used  [Ref.  11].  A  significant  effort  may  be  necessary  to  prevent  information 
loss  in  the  discretization  and  still,  it  may  be  prone  to  errors.  It  is  very  difficult  “to 
construct  an  efficient  preconditioning  method  (which  modifies  the  condition  number 
of  the  matrix  A)  for  arbitrary  ‘purely  algebraic’  problems  [Ref.  9].” 

What  we  really  want  to  do  is  to  solve  purely  algebraic  problems  with  the  same 
speed,  efficiency  and  reliability  that  multigrid  provides  to  geometrically  structured 
problems.  We  require  a  method  that  “should  use  only  information  in  the  matrix  of 
the  system  and  as  little  extra  information  as  possible  [Ref.  11].”  In  other  words. 
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we  want  to  extend  all  of  the  benefits  of  multigrid  methods  to  a  new  type  of  solver, 
algebraic  multigrid. 
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III.  ALGEBRAIC  MULTIGRID 


Algebraic  multigrid  (AMG)  methods  are  a  required  extension  of  the  multigrid 
methodology.  They  were  first  introduced  in  the  mid  1980’s  by  A.  Brandt  [Ref.  12],  S. 
McCormick  and  J.  Ruge  [Ref.  13]  and  further  developed  by  J.  Ruge  and  K.  Stuben 
in  their  paper,  Algebraic  Multigrid  [Ref.  14].  The  methods  arose  from  the  need  to 
solve  purely  algebraic  systems  (matrix  equations)  of  the  form 

n 

Au  =  f  or  =  fi  ii  =  l,...,n)  (HI.l) 

with  the  “multigrid-style”  cycling  process  without  any  underlying  geometry.  Al¬ 
though  there  is  nothing  special  about  the  form  in  (III-l)  as  is,  it  is  the  form  required 
for  our  work  in  AMG.  In  their  paper,  Ruge  and  Stuben  show  the  usefulness  of  AMG 
when  applied  to  a  “model  class”  of  matrices  -  symmetric  M-matrices.  (Recall  that 
the  matrix  A  is  an  M-matrix  if  an  >  0  and  aij  <  0  for  i  ^  j.)  They  also  relaxed 
the  assumption  to  weak  diagonal  dominance  (|ai,|  >  kiil?  i  =  1, . . . ,n)  and 

showed  uniform  two-level  convergence  for  that  case  as  well.  In  addition,  they  imply 
that  AMG  may  be  attractive  as  a  “black  box”  solver  for  algebraically  posed  elliptic 
problems  as  well  as  for  certain  other  types  of  operators  that  generate  matrices  with 
similar  characteristics. 

Moreover,  for  particular  types  of  matrix  equations,  algebraic  multigrid  is  a 
robust  and  efficient  matrix  equation  solver.  It  is  designed  to  solve  these  types  of 
matrix  equations  (III.l)  using  the  same  principles  of  standard  multigrid  methods, 
without  requiring  the  need  of  an  underlying  geometry  of  the  continuous  problem. 
Application  of  AMG  to  many  problems  involving  systems  that  generate  symmetric 
M-matrices  is  already  shown  to  be  efficient.  We  discuss  below  some  other  types  of 
problems  to  which  AMG  can  successfully  be  applied. 
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A.  WHY  ALGEBRAIC  MULTIGRID? 

AMG,  like  multigrid,  uses  the  same  multi-scale  hierarchy  to  solve  the  problem 
Au  =  /.  At  each  level  during  the  solution  process,  the  system  is  “solved”  and  the 
error  corrected  until  the  process  terminates  at  the  lowest  level  (when  no  more  levels 
are  required)  and  then  the  system  is  solved  by  using  one  of  a  host  of  direct  methods 
such  as  Gaussian  elimination.  Thus  far,  there  seems  to  be  little  difference,  but  shortly, 
we  will  see  that  there  is. 

It  is  a  solution  method  that  is  ideal  for  solving  more  general  large,  sparse 
algebraic  equations  because,  unlike  MG,  it  is  not  dependent  on  particular  domains 
or  operators.  Although  AMG  can  solve  standard  elliptic  problems  on  uniform  rect¬ 
angular  grids,  for  such  problems,  it  is  a  less  efficient  solver  than  highly  specialized 
geometric  multigrid  solvers.  However,  once  the  setup  phase  (discussed  in  the  next 
section)  is  completed,  AMG  is  quite  competitive.  But  when  geometric  grids  become 
impossible  to  design  and  MG  cannot  be  tailored  to  fit  the  problem,  this  is  when  AMG 
should  be  used.  It  has  been  shown  in  various  papers  such  as  [Ref.  13,  14,  15]  to  have 
favorable  results  for  these  more  complicated  domains. 

Even  though  originally  designed  to  solve  large,  sparse  matrix  equations  that 
involved  M-matrices,  application  of  AMG  to  other  types  of  not-so-sparse  matrices  has 
increased  in  recent  times.  This  is  primarily  due  to  current  problems  of  interest  which 
are  on  unstructured  grids  and  are  extremely  large,  usually  on  the  order  of  millions  of 
equations  and  unknowns.  In  fact,  in  some  particular  cases,  algebraic  multigrid  has 
been  proven  to  be  the  fastest  method  available.  One  problem,  for  example,  is  the 
transistor  placement  problem  in  integrated  circuit  design.  This  particular  problem 
concentrates  on  layout  optimization  and  solves  a  system  that  determines  the  optimal 
location  for  hundreds  of  thousands  of  cells,  which  consist  of  hundreds  of  transistors, 
on  the  chip  surface.  Regler  and  Rude  show  that  AMG  is  a  competitive  alternative  to 
current  methods  in  layout  optimization. 
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B.  APPLICATIONS  OF  AMG 

Algebraic  multigrid  methods  are  useful  when  it  comes  to  solving  certain  sys¬ 
tems  that  either  contain  no  geometric  structure  or  that  produce  irregular  grids  tailored 
to  complex  domains  or  operators.  When  the  application  of  standard  multigrid  meth¬ 
ods  proves  to  be  impossible  or  entails  a  high  degree  of  difficulty  on  certain  problems, 
algebraic  multigrid  may  be  a  better  choice.  These  types  of  problems  are  described  in 
more  detail  in  [Ref.  14]  but  here,  they  are  summarized.  In  general,  AMG  should  be 
considered  over  geometric  multigrid 

(1)  when  the  problem  involves  complex  domains  so  that  any  discretization  we 
choose  is  still  too  fine  to  serve  as  the  coarsest  grid.  In  this  case,  the  work  required 
in  geometric  multigrid  to  determine  the  interpolation  schemes  would  be  prohibitive. 
AMG  is  clearly  the  choice  here  since  the  coarsening  process  is  automatic.  In  fact,  for 
this  case,  AMG  performs  quite  favorably  in  terms  of  operation  count  and  CPU  time 
[Ref.  15]. 

(2)  when  uniform  coarsening  is  not  allowed  at  all  on  the  finest  grid.  This 
happens,  for  example,  when  a  finite  element  discretization  is  applied  using  irregular 
triangulations.  Since  AMG  requires  no  geometric  structure  (i.e.,  no  uniform  grids), 
it  is  ideal. 

(3)  when  the  interpolation  operator  is  chosen  so  that  it  becomes  very  difficult  to 
find  a  relaxation  process  that  complements  the  operator  in  smoothing  the  error  enough 
to  allow  a  decent  coarse  grid  correction.  Sometimes,  it  is  impossible  to  determine  a 
relaxation  process  at  all.  AMG  is  not  affected  here  since  interpolation  is  defined  by 
the  entries  in  the  matrix,  and  the  weights  are  chosen  tot  reflect  the  relative  strength 
of  each  connection. 

(4)  when  a  problem  is  entirely  discrete,  especially  if  the  problem  contains  no 
geometric  structure  at  all.  Since  AMG  does  not  depend  on  an  underlying  continuous 
problem,  it  is  clearly  the  method  of  choice.  This  applies  to  problems  where  we  are 
only  given  the  system  matrix  A  and  the  right-hand  side  vector  /. 
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C.  THE  STEPS  OF  AMG 

Algebraic  multigrid  methods  are  designed  with  regular  (geometric)  multigrid 
methods  in  mind,  in  that  they  use  the  principles  of  multigrid  but  do  not  require  nor 
rely  on  the  geometry  of  the  particular  problem  being  solved.  Instead,  they  explicitly 
use  information  from  the  matrix  of  the  system,  i.e.,  the  coefficients  of  the  unknowns. 
Application  of  AMG  to  the  problem  Au  =  f  requires  a  two-part  process.  There  is 
an  initial  step  which  defines  the  intergrid  transfer  and  coarse-grid  operators  and,  in 
addition,  chooses  the  coarse  grid  itself.  This  first  step  of  the  process  is  called  the 
setup  phase  and  most  of  the  work  involved  in  writing  the  algorithm  takes  place  in 
this  portion  of  AMG.  “The  generality  of  AMG  must  be  paid  for  by  a  setup  phase 
that  may  take  80%  or  more  of  the  overall  time  [Ref.  15]”.  The  second  part  of  the 
AMG  process  is  called  the  solution  phase.  In  this  part,  regular  multigrid  cycling  is 
performed  on  the  components  of  the  matrix  until  a  predetermined  tolerance  is  met. 

1.  The  Setup  Phase 

As  defined  in  [Ref.  14],  a  brief  outline  of  the  algorithm  used  in  the  setup  phase 
on  the  problem  Au  =  /  is  provided  below: 

1.  Set  m  =  1. 

2.  Choose  the  coarse  grid  and  define  interpolation  operator. 

3.  Set  and 

4.  When  is  small  enough,  set  ^  =  m  +  1  and  stop.  Otherwise,  let 

m  =  m  -h  1  and  return  to  step  2. 

Here,  we  define  q  to  be  the  number  that  represents  the  coarsest  grid.  This  number  is 
used  in  the  theoretical  bound  on  the  asymptotic  (V-cycle)  convergence  factor,  p.  In 
practice  though,  p  is  ^-independent  and  by  increasing  the  accuracy  of  interpolation 
for  smooth  errors,  we  find  that  AMG  is  quite  an  efficient  method.  Therefore,  the  large 
investment  in  writing  the  setup  phase  may  save  time  (and  cost)  in  the  setup  work  for 
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later  problems  that  use  the  same,  or  nearly  the  same,  matrix.  This  effectively  reduces 
the  time  required  for  the  solution  of  follow-on  problems. 

2.  The  Solution  Phase 

The  solution  phase  is  part  of  the  algorithm  in  which  regular  multigrid  cycling 
takes  place.  It  involves  a  relaxation  process  to  remove  the  oscillatory  modes  of  the 
error  and  an  intergrid  transfer  operators  to  move  the  problem  to  different  levels.  This 
phase  of  AMG,  when  the  operators  and  relaxation  are  chosen  properly,  is  a  very 
efficient  solver  for  the  problem  on  the  finest  grid  [Ref.  14]. 

D.  DEFINING  THE  “GRID” 

One  of  the  main  ideas  of  algebraic  multigrid  is  to  split  the  set  of  variables 
on  a  particular  level  into  two  disjoint  subsets.  The  first  set  contains  coarse  level 
variables  (C- variables)  and  the  second  is  the  complementaxy  set  of  fine  level  variables 
(F- variables).  The  fact  that  we  require  the  sets  to  be  disjoint  will  be  important  as  we 
will  see  shortly.  If  we  now  define  h  and  H  to  be  any  two  consecutive  levels  with  h  as 
the  fine  one  then  we  can  set  \  j  e  C}  and  =  {j  \  j  e  F}.  Furthermore, 

let  i  e  U  F^  be  defined  as  a  point  in  a  fictitious  plane.  In  this  sense,  i 

is  nothing  more  than  a  reference  to  the  variable  uf.  Now,  for  example, 
can  be  interpreted  as  a  grid  equation  on  a  fictitious  grid  0,^.  Furthermore,  we  define 
Qh  _  jph  qh  _  qH  (Jetail  in  defining  the  grid  can  be  quite  cumbersome 
and  it  is  not  our  intention  here  to  reproduce  what  can  be  found  in  [Ref.  14]. 

E.  CONNECTIONS  AND  CONVERGENCE 

Next  we  introduce  a  new  term  called  connections.  Connections  refer  to  the 
relationship  between  grid  points  in  the  sense  of  a  directed  graph  which  is  associated 
with  any  matrix.  On  any  level  h,  a  point  i  €  is  defined  to  be  (directly)  connected 
to  a  point  j  €  0^  if  Uij  ^  0.  Furthermore,  we  define  the  (direct)  neighborhood  of  a 
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point  i  by 


Nt  =  O'  €  OM  i  #  4  ^  0}-  (III.2) 

Let  us  now  concern  ourselves  with  interpolation  along  direct  connections.  In 
this  case,  we  consider  only  the  interpolation  points  Ci  with  (7,-  C  Ni  O  C  and  the 
corresponding  weights 

,  {ieF,  k£  Ci)  (III.3) 

with 

0  <  7/i  <  1/  ^  |4  .  (IIL4) 

leCi 

In  order  to  maintain  two-level  convergence,  it  is  sufficient  to  require  for  every  i  € 
F,k€  Ci 

0  <  a’y.wiu  <  13 14;,| ,  0  <  a^.(l  -  Si)  <  4 

i 

where  s,-  <  1,  which  is  ensured  by  (IIL4).  From  (III. 5),  we  can  easily  derive  more 
practical  conditions  to  effectively  develop  an  automatic  coarsening  algorithm  that  has 
^  as  an  input  parameter.  The  algorithm  will  also  choose  interpolation  with  weights 
(III.3)  satisfying  (IIL5). 

In  order  to  understand  the  effect  of  ^  in  (III.5),  we  state  the  following  result 
presented  in  [Ref.  14]. 

Theorem:  Let  yd  >  1  be  fixed.  Assume  for  any  symmetric,  weakly  diagonally 
dominant  M-matrix  the  C-points  are  picked  such  that,  for  each  i  G  F, 
there  is  a  non-empty  set  Ci  C.  NiC^C  with 

^  E  4  ^  4  (111.6) 

iiCi 

and  define  the  interpolation  weights  (III.3)  by  rji  =  l/Ej^Cy  o^ij-  Then  two-level 
convergence  is  satisfied. 

For  a  given  /?,  condition  (III.6)  is  satisfied  by  many  choices  of  the  sets  C 
and  Ci'i  but  regardless  of  how  we  chose  those  sets,  we  always  want  them  as  small  as 
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possible.  This  condition  requires  that  we  make  each  F-point  i  strongly  connected  to 
its  interpolation  points.  In  other  words,  we  actually  want  each  |ay|  [j  €  Ci)  to  be 
comparable  in  size  to  the  largest  of  the  la,*!  (fe  €  Ni)  or  equivalently, 

Icjjl  7max  |a,fc|  for  j  €  Ci  and  0  <  7  <  1. 

The  assumption  in  (III.6)  gets  weaker  as  ^  gets  larger;  but,  for  large  (III.6) 
allows  for  rapid  coarsening.  The  price  of  this  though  is  slower  two-level  convergence 
speed.  Although,  when  yd  =  1,  we  achieve  the  best  convergence  but  the  expense  of 
the  method  becomes  extreme.  Ruge  and  Stuben  explain  that  the  best  compromise 
is  when  yd  =  2.  This  is  when  about  50%  of  the  total  strength  of  connections  of  every 
F-point  will  be  represented  on  the  coarser  level. 
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IV.  CONSTRUCTING  THE 
INTERPOLATION  OPERATOR 


The  choosing  of  an  interpolation  operator  is  “dictated  by  the  desire  to  achieve 
good  theoretical  estimates  of  the  convergence  of  the  iterations  as  well  as  by  bounds 
on  the  computational  complexity  of  the  algorithm  [Ref.  11].”  It  is  of  the  essence 
in  defining  the  interpolation  operator  that  the  range  of  interpolation  approximates 
those  errors  not  effectively  reduced  by  relaxation.  To  make  it  successful  then,  AMG 
must  be  constructed  in  such  a  way  to  ensure  that  this  automatically  happens. 

A.  ALGEBRAIC  SMOOTHNESS 

For  the  problem  Au  =  /,  let  us  define  a  smoothing  process 
«new  =  G  U^U  +  I  -G  (A  )  / 

where  G  is  a  (linear)  smoothing  operator  such  that  G”*  :  T?.””*  7^”"* .  In  AMG,  an 

error  e  is  said  to  be  smooth  if  ||Ge||j  «  ||e||j.  In  other  words,  the  error  is  smooth  when 
it  is  slow  to  converge  with  respect  to  the  smoothing  operator.  If  we  now  assume  that 
e  is  a  smooth  error  then  the  residual  r  =  Ae  is  small  compared  to  e.  This  occurs  in 
most  of  the  common  relaxation  processes.  Knowing  this,  we  expect  |r,|  <C  an  |e,|  for 
all  i  and  so  for  (algebraically)  smooth  error  we  have  Ae  0.  A  good  approximation 
for  each  e,-  can  then  be  obtained  from 

Ae  =  anCi  +  ^  Cijej  ^  0 

j€Ni 

or 

e,-  =  —  X)  (IV- 1) 

da  jgjVi 

where  the  neighborhood  iV,-  of  point  i  is  defined  in  (III.2). 

With  this  fact,  it  is  easy  to  construct  an  interpolation  operator  which  guaran¬ 
tees  that  the  smooth  error  lies  in  the  range  of  interpolation.  One  steindard  operator 
might  be 
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ef  for  i  £  C  = 
for  i  e  F  = 

with  Wik  =  Sik  if  i  E  C'^.  Recall  that 


^ik  = 


for  i  ^  k 
for  i  =  k. 


However,  this  basic  “model”  interpolation  operator  isn’t  very  elFective  since  it 
is  not  usually  “local”  enough.  Here  we  define  “local”  to  be  Wik  =  0  for  all  i  € 
unless  k  E  C^,  where  Cf  C  C^.  For  reasons  of  efficiency,  we  require  that  Cf  be 
some  small  set  of  interpolation  points.  When  working  with  symmetric  M-matrices, 
[Ref.  14]  shows  for  the  case  when  |aij|  an  that  smooth  error  generally  varies 
slowly  in  the  direction  of  strong  connections.  This  fact  will  be  important  later  on 
when  constructing  weight  definitions  Wik  for  the  interpolation  operator  along  those 
connections. 


B.  INTERPOLATION  ALONG  DIRECT  CONNECTIONS 


Let  us  consider  only  those  interpolation  points  Ci  such  that  Ci  C  Ni  fl  C. 
The  interpolation  weights  then  have  the  form  Wik  =  rji  for  i  €  F  and  k  €  Ci 
where  0  <  77,-  <  |4|)~  .  We  require  this  definition  to  ensure  that  Wik  <  1. 

Let  us  take  a  look  at  the  case  when  the  matrix  A  is  strictly  an  M-matrix  (with  no 
assumption  on  symmetry),  i.e.,  an  >  0  and  aij  <  0  for  i  ^  j.  If  we  set  Ci  =  Ni,  then 
(IV.l)  gives 


e. 


(IV.2) 


Now  we  choose  Wik  =  \a'>k\  /an.  However,  this  would  imply  that  C  =  0,  but  we  want 
Ci  to  be  small;  so  generally,  we  desire  the  set  Di  =  Ni  —  Ci^^.  What  we  want  to  do 
is  to  “distribute”  the  non-interpolatory  connections  a^  (j  6  D/)  to  the  interpolatory 
point  an. 
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If  |aij|  is  small,  then  adding  a,j,  for  j  €  Di,  to  the  diagonal  element  an  should 
not  hurt  accuracy  very  much  since  the  contribution  from  Cij  is  negligible.  If,  however, 
\aij  I  is  large,  then  we  know  this  replacement  is  reasonable  since  smooth  error  generally 
varies  slowly  in  the  direction  of  strong  connections.  To  determine  the  error  in  this 
case,  we  begin  with  (IV.  1). 


then 


or 


an^i  —  ^  ^  ^ij^j 

jeNi 

=  —  ~ 

j€Ci  jeDi 

' - V - ' 

set  cj^ei 

~  dijCi 

jeCi  j€Di 


I  “I"  ^  ^ij^j 

\  jeDi  )  jeCi 

_ jeCi 

a>ii  "b  2^  ^ij 
jeDi 

_  jeCi _ 

^ii  "b  S 

jeDi 


(IV.3) 


Since  A  is  assumed  to  be  an  M-matrix,  the  denominator  in  (IV.3)  is  positive.  We 
now  have  a  good  representation  of  the  error. 


C.  THE  ITERATIVE  WEIGHT  DEFINITION 

Let  us  now  consider  the  entries  in  an  arbitrary  matrix  A.  Figure  12  shows  a 
graphical  relationship  between  various  points  in  the  matrix.  We  will  use  this  illus¬ 
tration  in  our  derivation  of  the  weights.  Recall  our  earlier  assumption  that  Ae  «  0 
when  e  is  smooth.  Then,  for  a  smooth  error  and  for  a  point  i  £  F,  we  have 

an^i  ^  ^  ^ik^k  ^  ^  a^jCj  ^  ^  (^V -4) 

k€Ci  jeF  leNi,  l£C,  l^Ci 
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l€  N; 


Figure  12.  Relationship  of  fine  grid  to  coarse  grid  points  in  the  matrix.  Points  i  and 
j  are  the  fine  grid  points;  the  others  are  coarse  grid  points,  both  weak  and  strongly 
connected,  interpolatory  and  non-interpolatory. 

and  similarly  for  a  point  j  £  F,  we  have 

keCj  ieF  meNj,  mec,  m^Cj 

To  make  the  interpolation  operator  “small,”  we  must  distribute  the  non-interpolatory 
connections  (aij  for  j  ^  Cf)  to  the  interpolatory  points.  In  basic  terms,  we  need  to 
approximate  ej  {j  €  F)  and  ei  (/  €  AT,-,  I  £  C,  I  ^  Cj)  in  terms  of  e*  (k  £  C,). 
Likewise, we  need  to  approximate  Cj  (i  £  F)  and  e^  (m  £  Nj,  m  £  C,  m  ^  Cj)  in 
terms  of  e*  (k  £  Cj). 

Let  us  assume  that  the  interpolation  weights  have  been  determined  for  the 
interpolatory  points  k  £  Cj.  That  is,  the  interpolatory  weights  to  define  Cj  for 
j  £  F  are  known.  Then  ej  could  be  approximated  by  a  weighted  average  (for  use  in 
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approximating  ej): 


so  in  general 


^jki^ki  ^jk2^k2 
Wjki  +  Wjk2 


^jk  ^k 

kSCi 

E  Wjk 

k&Ci 


(IV.6) 


Then  for  i  e  F,vre  substitute  ey  in  (IV.6)  into  ey  in  (IV.4).  This  leads  to 


Ctii&i  —  y  0,ik^k  ^  ^  Oj'y 
keCi  jeF 


i  E  yjjk 

\  keCi 


leNi,  lec,  liCi 


(IV.7) 


while  for  j  G  F,  we  replace  ti  in  (IV.5)  by  the  equivalent  approximation  for  Ci  in 
(IV.6).  This  gives  us 


^jj^j 


^jk^k  <^ji 

keCj  ieF 


i  E  ^ik&k\ 


keCj 


.  E  Wik 
\ 


) 


E 

mSNj,  m£C,  m^Cj 


(IV.8) 


If  we  look  at  the  two  previous  equations  closely,  we  see  that  (IV.7)  is  being 
used  to  define  Wik  which  are  used  in  (IV.8)  and  (IV.8)  is  being  used  to  define  Wjk 
which  are  used  in  (IV.7).  Thus,  the  interpolation  weights  at  point  i  are  dependent 
on  the  interpolation  weights  at  point  j,  and  vice-versa.  Hence,  we  have  an  implicit 
system  to  define  the  interpolatory  weights.  This  implicit  system  governs  how  the 
weights  are  computed  and  is  the  basis  for  the  iterative  weight  definition  (IWD)  of 
AMG. 


D.  INTERPOLATION  CONSTRUCTION  USING  THE 
ITERATIVE  WEIGHT  DEFINITION 

This  section  presents  an  algorithm  for  coding  the  iterative  weight  definition  of 
algebraic  multigrid.  We  begin  by  starting  at  a  fine  grid  point,  say  point  i,  and  dividing 
its  neighborhood,  Nj,  into  four  distinct  groups.  Figure  13  represents  an  illustration 
of  a  typical  neighborhood  of  an  arbitrary  point  i.  It  also  shows  the  elements  in  each 
of  the  groups  to  which  we  will  refer  in  this  discussion  of  construction. 
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We  begin  the  algorithm  by  dividing  the  entries  of  the  matrix  into  four  groups. 
For  each  point  i,  divide  Ni  into  four  groups: 

l€  N; 


Figure  13.  Illustration  of  the  neighborhood  (Ni)  of  a  fine  grid  point  i  and  its  relation- 


ship  to  other  fine  and  coarse-grid  points  in  the  matrix.  For  each  point  i,  Ni  is  divided 

into  four  distinct 

groups. 

Group  A 

Coarse-grid  interpolatory  points,  C,-.  In  Figure  13,  these  corres¬ 
pond  to  points  ^1,  and  k^. 

Group  B 

Coarse-grid  and  fine-grid  non-interpolatory  points  that  are  weak¬ 
ly  connected  to  point  i.  In  Figure  13,  these  are  points  m  and  /. 

Group  C 

Fine-grid  points  strongly  connected  to  point  i.  In  Figure  13,  this 

is  point  j. 

Group  D 

Coarse-grid  points  not  weakly  connected  to  point  i  but  not  in 

Ci.  In  Figure  13,  this  is  point  n. 
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1.  Initialization 

Step  1:  Initialize  the  variables  (vectors): 
For  all  i  €  F  do 
set 

Ci  =  Si  DC 


set 


define 


end 


Wik  = 


l€Ci 


=  {j  I  o,ij  is  small  } 


Here,  we  note  that  Si  is  the  set  of  points  strongly  connected  to  point  i  and  is  the 
set  of  points  weakly  connected  to  i. 

2.  Calculation 

Step  2:  Calculate  the  interpolatory  weights,  Wik,  by  using  the  approximation 
Ae  PH  0: 

For  each  e  €  F  do 
and  !=s 


'y  ]  ^ik^k 

(IV.9) 

k&Ci 

(IV.IO) 

(IV.ll) 

j€F,  j^Qi 

(IV.12) 

TlSOy  Tl^Ciy 

approximate  the  errors: 

for  group  B  points,  distribute  the  quantity  in  (IV.IO)  to  the 
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diagonal  by  the  approximation 


ei  «  Cj. 

for  group  C  points,  approximate 


(IV.13) 


E  WikCk 

keCi 

E  mk 

keCi 


(IV.14) 


but  distribute  the  quantity  in  (IV.  11)  to  the  diagonal  only 
when  Wik  >  0. 

for  group  D  points,  approximate 


E  OikCk 

keCj 

E  aik 

keCi 


(IV.15) 


but  distribute  quantity  in  (IV.12)  to  the  diagonal  only  when 
Oij  is  of  the  right  sign. 

Either  stop  here,  or 
For  each  i  E  F  do 
set 


Ci  =  {j  I  Wij  >  0} 

goto  Step  2. 

end 

end 

We  give  some  insight  to  the  iterative  nature  of  this  last  step.  The  routine  ends 
when  there  are  no  more  fine-grid  points  or  the  set  of  interpolatory  coarse-grid  points, 
Ci,  is  empty.  If  this  is  not  the  case,  we  check  the  next  fine-grid  point  i  and  refine  the 
set  Ci  to  include  only  those  weights  that  are  of  the  correct  (positive)  sign.  If  all  the 
weights  are  now  of  the  correct  sign,  then  we  end,  if  not  we  refine  Ci  again  and  repeat 
the  loop  until  we  meet  an  end  of  loop  criteria. 
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E.  INTERPRETATION  OF  THE  ALGORITHM 


In  this  section,  we  interpret  the  code  of  Section  D.  But  before  we  begin,  we 
make  a  few  simplifications.  Define  the  sets 

A  =  {kikeCi} 

B  =  {ilieQi} 

c  =  {jlJeFjia,} 

D  =  {n\n  €:C,n  ^  Ci,n  ^  flj} . 

Now,  the  equation  involving  (IV.9)-(IV.12)  becomes, 

a«ei  =  -  I  OikCk  +  ^  and  +  X]  H  )  .  (rV.16) 

\  A  B  C  D  / 

In  the  code,  the  sum  over  the  set  B  becomes  after  substitution  of  (IV.  13)  into  (IV.  10) 


'Eaud,  (IV.17) 

B 

the  sum  over  the  set  C  becomes  after  substitution  of  (IV.14)  into  (IV.ll)  assuming 

Wik  >  0 


(T.Wikek 


y  E  '^ik  '  ’ 

\keCi 


(IV.  18) 


and  the  sum  over  the  set  D  becomes  after  substitution  of  (IV.15)  into  (rV.12)  assum¬ 
ing  aij  is  of  the  right  sign 


D 


(  ^ik^k 
A _ 

1  ^ik 
\  keCi 


After  interchanging  sums  in  (IV.  18),  we  get 


E 


72  aij  \ 

C 


72  Wik 

Kk^Ci 


Wiktk 


(IV.19) 


(IV.20) 
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and  after  switching  the  sums  in  (IV.19),  we  get 


/  ^ 
E'  ^ 


I  ^ik 

\keCi 


^ik^k* 


(IV.21) 


With  (IV.17),  (IV.20)  and  (rV.21),  (IV.16)  becomes 


CliiCi  “f-  ^  ^  ^il ~ 
B 

or  equivalently 


/  XT 

E  ^ikek  +  E  ^ 

A  A 


(  E«f 


I  E  Wik 
\fc€C, 


WikCk 


A  ^ik 

\keCi 


dik^k 


Eoii  '( 

f  EOin 

■ 

«-E 

A 

aik  + 

c 

Wik  + 

D 

dik 

^  E  Wiki 

1  E  dik 

- 

\k€Ci  ) 

\keCi  / 

ek.  (IV.22) 


We  now  rewrite  (IV.22)  as 


ei  «E 

A 


X2®0 


<^ik  +  I  1  Wik  +  *  ^ 


JZ 


AGCj  > 


a*jb 


"h  ^ 

The  term  in  the  braces  is  Wik-  So  our  formal  definition  of  the  iterative  weight  definition 


f  ^k’ 


IS 
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V 


NUMERICAL  RESULTS 


In  this  chapter  we  present  some  numerical  results  on  problems  involving  dis¬ 
cretized  partial  differential  equations  as  well  as  non-geometrically  generated  matrix 
equations  using  algebraic  multigrid  methods.  As  will  be  seen  shortly,  AMG  methods 
using  the  iterative  weight  definition  (IWD)  may  not  be  robust  but  are  more  efficient 
on  certain  problems  than  standard  AMG  weight  definitions. 

Before  we  begin  with  the  results,  a  few  definitions  are  in  order.  We  define  a 
strong  connection  as 


\aij\  =  7 max  |aiA:|  for  j  e  Ci 

k€Ni 

where  7  is  a  parameter  that  varies:  0  <  7  <  1.  The  choice  we  use  in  our  experiments 
to  define  a  strong  connection  is  7  =  0.25  which  has  been  experimentally  shown  to 
yield  good  results.  Most  problems  presented  use  this  choice  of  7;  however,  for  some 
problems  we  test  with  7  =  0.5  because,  in  certain  cases,  favorable  results  have  been 
found  (and  is  another  avenue  of  research).  For  those  problems,  it  will  be  specifically 
stated  what  7  is.  For  reference  purposes,  (1,1)  V-cycles  with  Gauss-Seidel  relaxation 
and  C/F-ordering  of  points  were  used  in  all  tests.  The  convergence  factor  (p)  was 
computed  by 


where 

initial 


=  (IM\ 

vm; 


(V.l) 


n  equals  the  number  of  cycles  performed,  Rf  is  the  final  residual  and  Ri  is  the 
residual.  The  following  notation  will  also  be  used  throughout  this  chapter: 


cr^  Operator  complexity 
cr^  Grid  complexity 

p  the  asymptotic  convergence  factor  of  the  20‘^  V-cycle 

where  is  the  ratio  of  the  total  number  of  nonzero  entries  in  all  the  matrices  to  that 
of  the  number  of  nonzero  entries  in  the  fine-grid  matrix  and  cr^  is  the  ratio  of  the 
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total  number  of  points  on  all  grids  to  that  of  the  number  of  points  on  the  fine-grid, 
cr^  and  <7^  are  other  measures  that  describe  the  performance  of  AMG.  After  the  first 
few  experiments  we  use  p  instead  of  p  since  the  20*^  V-cycle  is  more  representative 
of  the  true  convergence  rate. 


A.  ISOTROPIC  PROBLEMS 

We  wish  to  establish  some  baseline  behaviors  of  each  method  in  the  absence 
of  any  irregularities.  We  discretize  the  Laplace  operator  (V^)  on  the  unit  square  with 
Dirichlet  boundary  conditions  using  several  different  finite  difference  discretizations. 
For  each  problem,  we  use  a  uniform  grid  with  mesh  size  h  =  1/64.  If  we  let  N  =  If  h 
then  we  create  a  square  matrix  of  size  x  (or  4096  x  4096  for  our  problems). 
The  stencils  we  use  are  as  follows: 


(1) 


A2 


-1 

-1  4  -1 

-1 


(3) 


1 


-1  -1  -1 

-1  8  -1 

-1  -1  -1 


(2) 

(4) 


1 


1 

20^2 


-1  -1  1 


1 

-ij 

-1 

-4 

-1 

-4 

20  -4 

-1 

-4 

-1 

These  stencils  arise  naturally  from  the  discretization  of 

-V‘^u{x,y)  =  f{x,y)  in  fl 
u  =  0  on  dfl. 

To  get  stencil  (1),  for  example,  we  replace  the  derivatives  of 


(V.2) 


Uxx  Uyy  —  (^-3) 

by  second  order  finite  differences.  If  we  let  =  1/N  and  hy  =  1/M,  where  N  and 
M  are  positive  integers  then 

'^xx  ^  ~  2uij  +  Ui+ij)  (V.4) 
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and 


‘^yy  ^  ^2  d"  •  (^•^) 

Substituting  (V.4)  and  (V.5)  back  into  (V.3)  and  with  the  definition  /ij  =  f{xi,yj), 
then  (V.2)  becomes 


^i-lj  d~  ‘^i+l,j 

hi 


1  d~  Uij^i 


=  /, 


Ui^j  —  '^NJ  —  'WijO  —  '^iyM  —  0 


(V.6) 


where  1  <  i  <  —  1  and  1  <  i  <  M  —  1. 

In  terms  of  a  stencil  then,  we  get 


-1 

“ 

-1 

2  2 

-1 

hi 

hi 

If  we  now  let  h  =  =  hy,  then  (V.7)  becomes 


-V2  = 


1 

h2 


-1 

-1  4  -1 

-1 


(V.7) 


which  is  exactly  what  we  were  trying  to  show.  The  other  stencils  can  be  similarly 
derived. 

Table  I  compares  the  results  of  the  iterative  weight  definition  against  standard 
AMG  weights  on  Laplace’s  equation  (I.l)  using  the  stencils  above.  In  all  tests  the  ini¬ 
tial  guess  for  t/(x,  y)  was  randomly  generated.  On  Laplace’s  equation  using  standard 
stencils  discretized  using  finite  differences,  we  note  that  there  is  little,  if  any,  differ¬ 
ence  in  p  of  the  iterative  method  over  the  standard  definition.  As  expected  though, 
AMG  produces  excellent  results  on  these  model  problems  using  either  method. 
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Table  I.  Convergence  rates  in  solving  Laplace ’s  equation  using  different  isotropic  op¬ 
erators  discretized  using  the  finite  difference  method. 


Standard  weights 

Iterative  weights 

Stencil 

P 

Complexity 

P 

Complexity 

(7^ 

'fmam 

1 

0.05037 

2.3714 

1.7009 

0.05037 

2.3726 

1.6990 

2 

2.2475 

1.6799 

0.07903 

2.2588 

1.6797 

3 

1.4162 

1.3464 

0.08109 

1.4172 

1.3440 

4 

0.06541 

3.7120 

1.9963 

0.07938 

3.9734 

2.0522 

B.  THE  ANISOTROPIC  PROBLEM  USING  THE  FI¬ 
NITE  DIFFERENCE  METHOD 

Now  that  we  have  some  baseline  results  of  both  methods  in  the  absence  of  any 
irregularities,  we  wish  to  try  a  new  problem, 

—euxx  —  Uyy  =  0  in 

(V.8) 

u  =  0  on  dU 

where  the  operator,  —eu^x  —  Uyy,  is  the  anisotropic  Laplacian.  Results  in  this  section 
are  given  for  various  levels  of  anisotropy  (e)  and  again  we  use  second  order  finite 
differences  to  discretize  the  problem.  This  section  will  demonstrate  the  ability  of 
AMG  to  tailor  the  coarsening  algorithm  to  the  individual  problem.  Other  results  can 
be  found  in  [Ref.  9].  The  domain  in  the  following  problem  is  again  the  unit  square 
with  Dirichlet  boundary  conditions.  The  problem  is  discretized  on  a  uniform  grid 
with  h  —  1/64,  using  the  5-point  stencil 

(V.9) 

We  will  now  apply  both  weight  methods  to  problem  (V.8)  with  /(a:,  y)  =  0. 
The  motivation  behind  this  test  is  to  have  knowledge  of  the  error.  With  the  error 


/l2 


-1 

—e  2-h2s  —e 
-1 
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known  exactly,  it  is  easier  to  compare  the  true  performance  of  both  methods.  As 
before,  the  problem  is  discretized  on  a  uniform  grid  with  h  =  l /64,  using  the  5-point 
stencil  in  (V.9)  above.  This  generates  a  4096  x  4096  matrix.  The  convergence  rate  in 
Table  II  is  that  of  the  20th  iteration  (p)  vice  the  average  defined  in  (V.l).  In  Table 
II,  we  see  that  the  iterative  weight  definition  produces  a  better  convergence  rate  over 
standard  AMG  weight  methods  for  all  values  of  e  in  the  anisotropic  problem  but  these 
improvements  are  utterly  insignificant  since  they  are  so  small.  In  addition,  IWD  does 
not  always  produce  better  results  in  terms  of  solver  complexity  or  grid  complexity. 
Clearly  the  standard  weight  definition  is  the  better  in  this  case  since  the  increased 
cost  associated  with  IWD  makes  it  unfavorable  when  the  results  are  practically  the 
same. 

Table  II.  The  anisotropic  problem  of  the  Laplacian  operator  discretized  using  the  finite 
difference  method. 


Standard  weights 

Iterative  weights 

€ 

P 

Complexity 

P 

2.6319 

1.9644 

1.9609 

0.0100 

0.06731 

2.7707 

1.9590 

miiki 

1.9590 

mm 

0.5000 

0.06349 

0.06085 

2.3272 

2.3714 

ESIiUsI 

1.6990 

2.4271 

1.7129 

1.6956 

3.5391 

1.9033 

0.05575 

3.5691 

1.9070 

100.00 

0.06753 

2.8047 

1.9658 

0.06720 

2.8043 

1.9653 

1000.0 

0.04544 

2.6319 

1.9644 

0.03946 

2.6299 

1.9609 

Again,  one  can  see  in  Table  II  that  the  differences,  if  any,  between  the  methods  are 
minimal.  For  these  simple  problems  presented  thus  far,  we  see  little  benefit  in  using 
the  iterative  weight  definition.  In  fact,  if  we  consider  the  cost  associated  with  IWD, 
we  conclude  that  the  method  is  worse  since  the  standard  weight  method  is  a  less 
expensive  routine. 
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C.  THE  ANISOTROPIC  PROBLEM  USING  THE  FI¬ 
NITE  ELEMENT  METHOD 

We  present  results  in  this  section  using  the  discretized  Laplacian  operator. 
Whereas  Sections  A  and  B  involved  discretizations  using  the  finite  difference  method, 
we  apply  the  finite  element  method  to  see  if  there  are  any  improvements. 

Table  III  provides  results  generated  by  elongating  the  discretized  grid  for 
Laplacian  operator  in  the  a:-direction.  We  also  compare  the  convergence  rates  of 
both  methods  as  we  change  the  definition  of  strong  connection  from  7  =  0.25  to 
7  =  0.5.  In  the  table,  dx  =  edy  means  that  for  every  point  in  the  y-direction  (dj/), 
there  corresponds  e  more  points  in  the  a:-direction  {dx).  The  intention  of  this  test 
is  to  see  how  well  the  iterative  method  tailors  to  semi- coarsening  to  improve  the 
convergence  rate. 


Table  III.  The  anisotropic  problem  with  zero  RHS  discretized  using  the  finite  element 
method. 


standard  weight 

iterative  weight 

size  {N  X  N) 

dx  =  edy 

7  =  0.25 

7  =  0.5 

7  =  0.25  7  =  0.5 

matrix 

grid 

e 

P 

400 

20 

25 

0.14 

0.14 

0.14 

0.14 

900 

30 

10 

0.45 

0.23 

0.13 

0.13 

900 

30 

100 

0.83 

0.53 

0.82 

0.23 

1225 

35 

50 

0.25 

0.14 

0.15 

0.14 

10000 

100 

10 

0.47 

0.24 

0.14 

0.14 

10000 

100 

100 

0.93 

0.55 

0.93 

0.28 

From  inspection  of  Table  III,  it  is  clear  that  for  certain  problems,  the  iterative 
weight  definition  is  significantly  better.  In  some  cases,  IWD  is  better  by  over  |  and 
the  extra  cost  associated  with  using  this  method  has  more  than  paid  for  itself  with  the 
large  reduction  in  convergence  rate.  The  natural  questions  then  are  why  is  it  better 
on  some  and  not  on  the  others?  What  are  the  special  properties  of  the  problems  in 
which  IWD  outperforms  the  standard  weight  definition?  Does  semi-coarsening  play 
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a  role?  Before  we  can  answer  these  questions,  we  look  at  some  other  problems  with 
domains  that  are  much  more  complicated. 

D.  PROBLEMS  WITH  COMPLEX  DOMAINS 

In  this  section,  we  present  some  results  for  some  complicated  domain  problems 
using  the  standard  Laplacian  operator 

d  d 

dx^  dx^‘ 

As  in  Sections  C,  we  use  the  finite  element  method  to  discretize  the  problem.  Figures 
14,  15  16,  17,  18  and  19  are  graphic  illustrations  (meshes)  of  the  problems  we  inves¬ 
tigate.  The  size  of  each  problem  varies  and  is  noted  in  the  following  tables.  Table 
IV  shows  the  results  of  some  initial  experiments  on  these  more  complicated  domains. 
Much  like  the  problems  we  have  already  seen,  we  expect  that  little  will  be  gained  by 
using  the  iterative  weight  definition  since  on  standard  Laplacian  operators,  regular 
AMG  already  works  great. 

1.  The  Meshes 


Figure  14.  Mesh  1. 


51 


Figure  15.  Mesh  2. 


Figure  17.  Mesh  4- 


Figure  18.  Mesh  5. 
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Figure  19.  Mesh  6. 

2.  Results  of  Complex  Domain  Problems 


Table  IV.  Convergence  rates  of  problems  high  in  complexity  using  the  standard  Lapla- 
cian  operator  discretized  using  the  finite  element  method. 


Standard  weights 

Iterative  weights 

Mesh 

Size 

P 

Complexity 

P 

Complexity 

(7^ 

<7^ 

1 

561 

0.8717 

2.1990 

1.7184 

0.8496 

2.1998 

1.7184 

2 

2.1250 

1.6530 

0.0902 

2.1118 

3 

321 

0.0746 

3.0207 

■liMfci 

4 

757 

0.1386 

2.9944 

1.9511 

0.1779 

2.9780 

1.9406 

1.9936 

0.1460 

3.3000 

1.9745 

6 

1080 

0.1501 

2.3064 

1.7630 

0.1399 

2.3541 

1.7769 

With  these  results,  we  suppose  that  with  simple  operators  such  as  the  standard 
Laplacian  or  even  the  slightly  skewed  Laplacian  the  iterative  weight  definition  will 
not  perform  better  than  regular  AMG  methods.  This  is  largely  due  to  the  fact  that 
conventional  MG  works  extremely  well  on  standard  (elliptic)  operators.  We  must  test 
our  new  method  using  a  more  irregular  type  of  operator  to  see  if  IWD  will  fare  better 
than  standard  weight  definitions. 
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E.  COMPLEX  DOMAINS  AND  OPERATORS 


Laplacian  operators  have  not  seemed  to  produce  results  in  favor  of  the  iterative 
method  so  we  hypothesize  that  IWD  may  outperform  the  standaxd  weight  definition 
method  when  the  operator  is  more  complex.  To  get  a  better  idea  of  what  types  of 
operators  or  on  what  types  of  problems  the  iterative  weight  definition  improves  p, 
we  develop  a  more  complex  operator  and  reevaluate  some  of  the  same  problems  as  in 
Section  D. 

Specifically,  we  experiment  with  those  geometries  as  illustrated  in  Figures  14, 
15,  17,  and  16.  The  size  of  each  problem  here  also  varies  and  may  be  different  from 
those  presented  in  Table  IV.  The  sizes  of  each  problem  are  included  in  the  tables. 
For  problems  in  this  section,  we  use  the  irregular  operator 


Table  V  provides  results  of  both  methods  using  the  definition  of  strong  con¬ 
nection  7  =  0.25  whereas  Table  VI  uses  7  =  0.5  as  the  definition. 


Table  V.  Convergence  rates  of  complex  problems  using  an  irregular  operator  discretized 
using  the  finite  element  method  with  7  =  0.25. 


Standard  weights 

Iterative  weights 

Mesh 

Size 

A 

P 

Complexity 

A 

P 

Complexity 

1 

561 

0.1367 

2.3281 

1.7558 

0.1432 

2.3891 

1.7825 

513 

0.0984 

1.9415 

0.1279 

2.3543 

1.9376 

3 

1249 

0.5912 

2.0168 

4 

2889 

0.7949 

0.8880 

3.2519 

1.9740 

Tables  V  and  VI  and  reveals  something  interesting.  The  iterative  weight  defi¬ 
nition  significantly  improves  the  convergence  rate  over  the  standard  weight  definition 
by  approximately  |,  but  only  in  one  case.  What  is  special  about  that  case?  Clearly, 
it  is  not  this  irregular  operator  alone  nor  the  specific  problem  but  a  combination  of 
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Table  VI.  Convergence  rates  of  complex  problems  using  an  irregular  operator  dis¬ 
cretized  using  the  finite  element  method  with  7  =  0.5. 


Standard  weights 

Iterative  weights 

Problem 

Size 

P 

Complexity 

P 

Complexity 

(7^ 

1 

561 

0.1194 

2.7501 

1.9305 

2 

513 

0.2863 

3 

1249 

0.7573 

BtilMil 

4 

2889 

0.9013 

2.9619 

2.0142 

0.9400 

2.9612 

2.0059 

both  that  holds  the  answers  to  our  quest  for  classification.  Moreover,  in  the  other 
three  problems,  IWD  performs  worse. 

Let  us  now  investigate  the  case  (mesh  3)  where  IWD  works  the  best  to  see  if 
there  are  any  clues  to  its  usage.  Figure  20  shows  not  only  the  sparsity  pattern  of  the 
matrix  generated  but  more  to  the  point,  the  locations  of  the  entries  in  the  matrix. 


Figure  20.  The  sparsity  pattern  of  mesh  3. 

Figure  21  shows  where  the  positive  and  negative  entries  are.  Together,  these 
figures  reveal  that  mesh  3  is  not  a  M-matrix  and  so  we  have  a  non-M-matrix  case  where 
the  iterative  weight  definition  should  be  used  over  regular  AMG  weight  definitions. 


Figure  21.  Location  of  positive  and  negative  entries  of  mesh  3. 


Table  VII  shows  the  difference  between  the  iterative  weight  definition  and  the 
standard  weight  method  over  a  sequence  of  refinements  on  that  particular  problem 
(mesh  3).  In  all  cases,  we  see  a  major  improvement  in  convergence  rate. 


Table  VII.  Convergence  rates  of  the  problem  in  Figure  16  using  an  irregular  operator 
discretized  using  the  finite  element  method.  Different  refinements  to  the  mesh  are 
presented. 


Standard  weights 

Iterative  weights 

Size 

P 

Complexity 

P 

Complexity 

1 

85 

0.3967 

2.6373 

1.9176 

0.3411 

2.7199 

1.9412 

2 

0.4368 

retiiEEi 

0.2684 

3.4382 

2.0467 

3 

1249 

0.5912 

3.6256 

QBI 

0.2877 

3.6133 

2.0168 

4 

4929 

2.0016 

0.3582 

3.6430 

1.9846 

Table  VIII  shows  other  experiment  involving  this  irregular  operator  on  the 
problem  (mesh  5)  shown  in  Figure  18. 
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Table  VIII.  Convergence  rates  of  the  problem  in  Figure  18  using  an  irregular  operator 
discretized  using  the  finite  element  method.  Different  refinements  to  the  mesh  are 
presented. 


Standard  weights 

Iterative  weights 

Refine 

Size 

P 

Complexity 

P 

Complexity 

330 

1.0000 

2.9523 

1.9940 

0.9906 

2.9927 

2.0121 

1 

1256 

0.9888 

3.3843 

2.0677 

3.2193 

2.0239 

2 

4896 

0.8459 

3.6010 

2.0666 

3.4574 

2.0345 

In  this  case,  both  methods  have  seemed  to  fail  all  together.  We  conclude  here  that 
the  operator  alone  isn’t  the  cause  for  the  improvements  seen  in  Table  VII  in  using 
IWD.  It  must  be  a  combination  of  both  the  geometry  and  the  operator. 

F.  CONCLUDING  REMARKS 

The  accuracy  and  efficiency  of  AMG  in  solving  systems  involving  symmetric 
M-raatrices  has  been  shown  in  many  papers.  The  use  of  AMG  for  solving  other  types 
of  matrix  equations  is  a  topic  of  extensive  research  and  variations  in  computational 
methodology  during  the  interpolation  phase  have  resulted  in  differing  convergence 
rates.  We  have  found  that  regular  AMG  interpolation  weight  definitions  are  inade¬ 
quate  for  solving  certain  discretized  systems  that  do  not  lead  to  M-matrices.  For  these 
types  of  problems,  an  iterative  approach  to  determining  the  interpolatory  weights  was 
useful. 

In  applying  the  iterative  interpolatory  weight  definition  of  AMG,  we  have  care¬ 
fully  balanced  the  requirement  of  keeping  the  interpolatory  set  of  points  “small”  in 
order  to  reduce  solver  complexity  while  at  the  same  time,  maintaining  accurate  inter¬ 
polation  to  correctly  represent  the  coarse-grid  function  on  the  fine  grid.  In  addition, 
the  extra  work  involved  in  using  IWD  does  not  significantly  add  to  setup  phase  costs. 

Experimental  results  have  shown  that  the  iterative  weight  definition  does  not 
significantly  improve  the  convergence  rate  over  standard  AMG  when  applied  to  M- 
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matrices,  which  we  anticipated.  However,  the  improvement  becomes  significant  when 
solving  certain  types  of  complicated,  non-standard  algebraic  equations  although  it  is 
unclear  at  this  stage  of  development  what  details  are  required  to  cause  the  iterative 
weight  definition  to  outperform  the  regular  weight  definitions. 

We  have  seen  that  IWD  is  not  robust;  however,  there  are  specific  problems  on 
which  IWD  should  be  used.  When  using  a  standard  operator  such  as  the  Laplacian, 
regular  AMG  should  be  used  since  its  performance  is  superb  and  it  is  less  expensive 
to  use  than  the  iterative  weight  definition  we  presented  here.  However,  when  the 
operator  becomes  irregular  and  the  domain  more  complex,  IWD  has  been  shown  to 
dramatically  improve  the  convergence  rate  over  current  AMG  weight  definitions  and 
should  be  considered  a  viable  option  in  the  future. 
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VI.  FUTURE  RESEARCH 


Although  we  have  shown  that  the  iterative  weight  definition  (IWD)  of  AMG 
isn’t  robust,  we  have  found  certain  problems  where  it  has  improved  the  convergence 
rate  significantly.  We  are  currently  searching  for  more  data  that  will  enable  us  to 
classify  the  types  of  problems  on  which  IWD  works  best.  There  are  still  many  avenues 
of  research  and  IWD  is  just  one  such  direction.  Algebraic  multigrid  research  is  open 
to  a  plethora  of  new  and  interesting  ideas.  The  need  for  more  efficient  and  robust 
solvers  is  never-ending  and  is  the  driving  force  that  brought  about  AMG  in  the  first 
place.  To  discuss  all  the  new  possibilities  of  improved  AMG  would  be  a  paper  in 
itself;  therefore,  we  restrict  our  attention  to  that  of  just  a  few  new  methods  on  the 
forefront  of  improved  interpolation. 

A.  INTERPOLATION  USING  EIGENVECTORS 

One  method  that  can  improve  interpolation  is  based  on  the  idea  of  eigenvec¬ 
tor  approximation  and  was  originally  presented  in  1985  by  Ruge  and  Stuben  in  [Ref. 
14].  Approximation  using  eigenvectors  is  a  method  used  to  modify  interpolation  and 
the  coarse  grid  matrices  on  all  levels.  The  basis  of  the  idea  is  that  the  eigenvectors 
corresponding  to  the  smallest  eigenvalues  are,  in  general,  the  algebraically  smoothest 
functions  and  so  must  be  better  approximated  by  the  range  of  interpolation.  Unfor¬ 
tunately,  the  eigenvector  approach  requires  the  computation  of  the  smallest  eigenval¬ 
ues  and  their  corresponding  (approximate)  eigenvectors.  However,  there  are  efficient 
methods  for  finding  them  and  so  it  may  not  be  as  computationally  expensive  as  one 
would  think.  Computational  accuracy  in  determining  them  is  raxely  needed  since 
linear  interpolation  (under  the  assumption  that  smooth  functions  axe  locally  linear) 
“usually  ensures  accurate  enough  interpolation  of  the  needed  eigenvectors  when  stan¬ 
dard  interpolation  does  not  [Ref.  14].” 
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To  make  the  method  efficient,  we  must  integrate  the  computation  of  the  eigen¬ 
vector  with  the  updating  of  the  interpolation.  Before  the  coarse-grid  correction  can  be 
applied  on  each  level,  we  must  update  the  interpolation  and  the  coarse-grid  operator 
with  the  current  eigenvector  approximation.  This  must  happen  before  the  operation 
can  continue  to  the  next  coarser  grid.  Ruge  and  Stuben  note  that  this  can  be  done 
efficiently  and  that  several  eigenvector  approximations  can  be  calculated  simultane¬ 
ously,  if  required.  The  concept  of  the  eigenvector  approximation  is  a  tool  used  only 
for  the  improvement  of  interpolation,  after  which  we  use  AMG  in  the  standard  way 
to  solve  the  problem. 

B.  THE  LEAST-SQUARES  IDEA 

Tom  Manteuffel  (University  of  Colorado  at  Boulder)  and  Van  Henson  (Naval 
Postgraduate  School)  recently  introduced  a  new  method  that  incorporates  a  least 
squares  solution  on  the  local  level.  Initial  trials  have  produced  quite  favorable  results. 
The  idea  behind  this  new  method  involves  an  singular  value  decomposition  (SVD)  of 
the  matrix  on  a  local  level.  Presented  here  is  a  brief  outline  of  how  it  works. 

1.  Let  i  be  a  fine  point.  Define  Ni  to  be  the  set  of  points  that  includes  the 
fine  point  i  and  those  points  p  connected  to  point  i.  Furthermore,  define  the  set  N2 
to  include  the  set  Ni  and  all  the  points  q  connected  those  points  in  Ni  and  continue 
this  process  (e.g.,  Nz,  N4, . . .)  until  there  are  no  points  remaining. 

2.  Now  select  from  the  matrix  A  the  rows  corresponding  to  A^i.  Remove  all 
columns  that  contain  only  zeros.  This  removal  of  zero-columns  yields  a  n  x  m  matrix 
which  we  call  S.  For  example,  if  we  start  with  a  nine-point  stencil  on  regular  grid, 
then  we  end  up  with  a  9  x  25  matrix,  S. 

3.  Compute  S  =  WEV,  the  singular-value  decomposition  of  S.  At  least  the 
last  m  -  n  (or  16,  using  our  nine-point  stencil)  columns  of  V  are  null-space  vectors  of 
S.  Let  X  be  the  m  x  (m  —  n)  matrix  (e.g.,  25  x  16)  comprising  these  m  —  n  columns 
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of  y. 


4.  Let  T  be  the  m  x  m  matrix  whose  rows  correspond  to  the  same  points  as 
do  the  columns  of  S.  In  the  example  of  the  nine-point  stencil,  T  is  of  size  25  x  25.  T 
is  the  matrix  associated  with  solving  the  homogeneous  Dirichlet  problem  on  N2. 

5.  We  now  reorthogonalize  the  columns  of  X  to  be  [T-orthogonal,  using  a 
Gram-Schmidt  process  with  respect  to  the  T-inner  product 

j-i 

=  {Txj,  Wk)  Wk  for  j  =  1,  2, . . . ,  m  -  n, 

k=i 

where 

Uj 

Uj) 

and  define  W  be  the  mx  (m  —  n)  matrix  whose  columns  axe  u?i,  W2,  tos,  . . . ,  Wm-n- 

6.  Select  I,  a  set  of  k  points  to  be  the  interpolatory  points  from  which  the 
value  at  i  is  to  be  interpolated.  Normally  (but  not  always)  the  points  in  I  are  chosen 
from  among  the  C-points  connected  to  i.  Permute  the  rows  of  W  so  that  the  first  k 
rows  contain  the  values  of  the  reorthogonalized  singular  vectors  corresponding  to  the 
points  in  I.  The  (^-1-1)**  row  contains  the  values  of  the  singular  vectors  corresponding 
to  the  point  i. 

7.  Perform  k  Householder  reflections  from  the  right  to  bring  the  matrix  (e.g., 
25  X  16)  W  into  the  form: 

*  0  0 

*  *  0  0 

*  *  *  0  0 

:|c  sf:  5|C  5jc  5jc  .  .  .  jf: 

*  •  •  •  * 

This  reflection  does  not  alter  the  span  of  the  columns,  since  the  Householder  reflection 
is  a  unitary  operation. 


*—  row  for  interpolation  point 
row  for  interpolation  point 
<—  row  for  interpolation  point 
•e—  row  coressponding  to  point  i 
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8.  Now  determine  a  linear  combination  of  the  first  k  rows  that  gives  the  first 
k  entries  in  the  (A:  +  1)®*  row.  The  weights  of  that  linear  combination  are  the  desired 
interpolation  weights. 

Note  that  the  size  of  the  entries  in  the  {k  +  1)**  row  beyond  the  k*^  column 
give  a  “residual”  measure  of  how  well  the  process  works.  If  all  the  entries  are  “small” 
then  the  weights  do  a  good  job  approximating  the  singular  vectors  at  point  i.  If  some 
of  them  are  “large”  then  we  may  need  to  add  another  interpolation  point  to  the  set 
or,  the  set  of  equations  selected  in  step  2  should  be  expanded  to  be  those  in  N2  vice 
those  in  Ni. 

Initial  experiments  with  this  method  have  produced  excellent  results  on  several 
standard  problems.  On  certain  problems  requiring  semi-coaxsening  (i.e.,  different 
coarse-grid  spacings  in  different  coordinate  directions),  this  is  the  only  method  known 
to  produce  the  “correct”  interpolation  weights.  However,  the  method  is  extremely 
expensive,  requiring  small  SVD  calculations  at  every  fine-grid  point.  Achieving  the 
robustness  of  this  approach  at  lower  cost  is  a  major  focus  of  ongoing  research. 

C.  THE  COMPOSITE  GRID  FORMULATION 

Steve  McCormick  (University  of  Colorado  at  Boulder)  has  presented  some  ini¬ 
tial  work  on  improving  interpolation  by  approximating  “globally  smooth”  components 
in  the  neighborhood  of  the  interpolation  region.  His  idea  is  to  use  composite  grids  to 
“solve”  the  local  problem  on  the  coarse  grids.  What  he  recommends  is  to  have  the 
target  point  i  and  all  its  neighbors  on  the  local  fine  grid,  twice  removed  neighbors  on 
the  next  coarser  and  more  global  grid,  and  so  on,  then  use  the  current  AMG  coars¬ 
ening  to  solve  the  local  problem  which  could  in  fact  just  be  the  least-squares  scheme 
described  in  Section  B.  In  the  following  discussion,  L  is  the  n  x  n  matrix  AMG  is 
given  to  “invert”  and  I  is  the  q  x  q  identity  matrix. 
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Let  us  assume  that  a  point  i  is  to  be  interpolated  from  q  C  neighbors.  Define  a 
composite  grid  consisting  of  the  fine-grid  in  the  interpolation  region  and  progressively 
coarser  grids  as  we  move  away  from  the  target  point  i.  We  now  use  AMG  on  the 
composite  grid  with  the  current  coarsening  to  compute  an  n  x  5  matrix  Q  with  the 
following  property: 

Q  minimizes  p(Q'^LQ) 

subject  to 

(PQfPQ  =  I 

where  p{A)  denotes  the  spectral  radius  of  A.  Note  that  the  qxq  matrix  PQ  denotes 
the  restriction  of  Q  to  the  q  interpolatory  C  points. 

The  motivation  here  is  to  find  the  vectors  that  lie  in  the  near  null  space  of  L 
that  are  very  distinct  on  the  C  points.  Included  in  the  space  determined  by  Q  (i.e., 
its  span)  is  the  “eigenvector”  p  of  L  belonging  to  the  smallest  eigenvalue  where  the 
local  normalization,  {Pp)^Pp  =  1,  is  used. 

This  scheme  is  reasonably  cheap  except  that  the  usual  coarse  grids  would 
introduce  O(log  n)  complexity.  McCormick  notes  that  all  of  the  coarse  grids  may  not 
be  needed  for  the  interpolatory  region  and  is  hopeful  that  we’ll  need  just  a  few  very 
coarse  ones.  In  this  manner,  we  can  get  the  interpolation  we  want  in  a  fairly  general 
setting. 

If  we  simplify  this  local  problem  by  interpreting  the  columns  of  Q  to  be  vectors 
truly  defined  on  the  fine-grid  then  it  should  be  fairly  easy  to  implement.  With  only 
one  target  point,  it’s  easy  to  see.  Then,  Q  is  an  n  vector.  If  we  now  partition  Q  into 
its  components  u,  belonging  to  the  C  points,  and  v  belonging  to  all  other  points,  then 
the  constraint  (PQ)^PQ  =  I  just  becomes  p^p  =  1. 

If  we  partition  L  corresponding  to  the  u-v  partition,  we  have 

A  B 
C 


65 


and  it’s  can  be  seen  that  the  local  problem  is  just  the  eigenvalue  problem: 


minimize  RQ{u) 

where 

Du 

RQ{u)  = 

u-'u 

and 

D  =  A- 

Note  that  this  determines  u,  and  then  v  can  be  computed  by  u  =  —C~^B^u.  Now, 
when  we  want  more  targets,  we  just  solve  for  more  eigenvectors  of  D.  In  fact,  we 
want  all  q  eigenvectors  of  the  q  x  q  matrix  D.  Once  these  eigenvectors  are  known,  a 
least-squares  procedure  could  be  used  to  select  interpolation  weights,  similar  to  the 
manner  described  in  Section  B.  Preliminary  experiments  with  the  method  outlined 
in  this  section  have  yielded  some  favorable  indications,  but  there  is  much  yet  to  be 
accomplished  before  the  efficacy  of  the  approach  can  be  demonstrated. 

The  three  avenues  of  research  outlined  here  are  all  active  areas  of  effort.  All 
are  focused  on  the  principle  that  accurate  interpolation  of  “local”  null-space  and 
near  null-space  vectors  will  ensure  accurate  approximation  of  the  global  equivalents. 
Moreover,  an  increase  in  accuracy  of  interpolation  for  smooth  errors  results  in  a  more 
robust  and  efficient  V-cycle  convergence  making  algebraic  multigrid  a  competitive 
alternative  to  current  methods.  In  fact,  in  some  cases,  there  is  no  better  matrix 
solver.  Even  though  this  research  is  in  its  infancy,  indications  now  point  toward 
eventual  success  in  using  these  new  ideas  to  produce  a  robust,  accurate,  and  efficient 
interpolation. 
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GLOSSARY 


algebraic  multigrid  a  numerical  method  used  to  solve  particular  algebraic  (linear) 
systems  with  the  “multigrid-style”  cycling  process 

boundary  conditions  a  condition  imposed  on  a  bounding  surface  (in  three  dimen¬ 
sions)  or  line  (in  two  dimensions)  or  at  a  bounding  point  (in  one  dimension) 
coarse-grid  operator 

connections  a  point  i  is  connected  (directly)  to  a  point  j  €  if  ^  0 

diagonal  dominance  (strict)  k«i  >  Mij I  i  — 

diagonal  dominance  (weak)  |aii|  >  Wij\  i  =  1, . . . ,  n 

Dirichlet  boundary  conditions  u  =  0  on  the  boundary 

injection  one  type  of  a  linear  restriction  operator 

iterative  methods  a.k.a  relaxation 

Fourier  modes  vj  =  sin  where  k  is  the  wave  number 

Galerkin  condition  known  as  the  coarse-grid  operator 

interpolation  operator 
Laplace’s  equation  =  0 

M-matrix  a  matrix  which  is  positive  definite  and  whose  off-diagonal  components 
are  nonpositive 

Poisson’s  equation  V^u  =  5,  5  7^  0 

restriction  operator  7?^””*+’ 

smoothing  operator  G  :  72^”  — >  72",  G  acts  on  e  to  remove  the  oscillatory  modes 
of  the  error 

smoothing  property  the  property  of  an  iterative  method  that  quickly  removes 
the  oscillatory  modes  but  leaves  the  smooth  (or  low-frequency)  modes  of  the  error 
strongly  connected  a  point  i  is  strongly  connected  to  a  point  j  if  \aij\  ~  max|ajfc| 
for  j  £  Ci  and  k  €  Ni 

variational  properties  7^+^  =  (7^^^)^  and  =  7^+M’"7^^i 
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wave  number  the  fcth  mode  consists  of  A;/2  full  sine  waves  on  the  interval  (0,7V) 
C- variables  coarse  level  variables 
F- variables  fine  level  variables 


68 


LIST  OF  REFERENCES 


[1]  Achi  Brandt.  Multi-level  adaptive  technique  (MEAT)  for  fast  numerical  solution 
to  boundary  value  problems.  In  Proceedings  of  the  3rd  International  Conference 
on  Numerical  Methods  in  Fluid  Mechanics,  volume  1,  pages  82-89,  1973. 

[2]  Pieter  Wesseling.  An  Introduction  to  Multigrid  Methods.  John  Wiley  and  Sons, 
New  York,  New  York,  1992. 

[3]  F.  Krawczyk  and  T.  Weiland.  Use  of  a  multigrid  solver  in  the  mafia  module 
S3  for  electro-  and  magnetostatic  problems.  IEEE  Transactions  on  Magnetics, 
26(2);747-750,  March  1990. 

[4]  Van  Henson,  et.  al.  Multilevel  image  reconstruction  with  natural  pixels.  In 
SIAM  Journal  on  Scientific  Computing,  volume  17,  pages  193-216,  Philadelphia, 
Pennsylvania,  January  1996.  SIAM. 

[5]  Tin-Su  Pan  and  Andrew  E.  Yagle.  Numerical  study  of  multigrid  implementa¬ 
tions  of  some  iterative  image  reconstruction  algorithms.  IEEE  Nuclear  Science 
Symposium  and  Medical  Imaging  Conference,  3:2033-2037,  November  1991. 

[6]  P.  Vanek,  J.  Mandel,  and  M.  Brezina.  Algebraic  multigrid  by  smoothed  aggrega¬ 
tion  for  second  and  fourth  order  eUiptic  problems.  In  N.  D.  Melson,  T.  A.  Man- 
teuffel,  S.  F.  McCormick,  and  C.  C.  Douglas,  editors.  Seventh  Copper  Mountain 
Conference  on  Multigrid  Methods,  volume  CP  3339,  pages  721-735,  Hampton, 
VA,  1996.  NASA. 

[7]  William  L.  Briggs.  A  Multigrid  Tutorial.  SIAM,  Philadelphia,  Pennsylvania, 
1987. 

[8]  Gene  Golub  and  Charles  Van  Loan.  Matrix  Computations.  The  John  Hopskins 
University  Press,  Baltimore,  Maryland,  1989. 

[9]  G.  Golubovici  and  C.  Popa.  Interpolation  and  related  coarsening  techniques  for 
the  algebraic  multigrid  method.  In  P.  W.  Hemker  and  P  Wesseling,  editors. 
Multigrid  Methods  IV,  number  116  in  International  Series  of  Numerical  Mathe¬ 
matics,  pages  201-213,  Basel,  Switzerland,  July  1994.  ISNM,  Birkhauser  Verlag. 

[10]  Achi  Brandt.  Multigrid  techniques:  1984  guide.  Computational  Fluid  Dynamics 
Lecture  Series,  March  1984. 

[11]  Petr  Vanek,  Jan  Mandel,  and  Marian  Brezina.  Algebraic  multigrid  on  un¬ 
structured  meshes.  Technical  Report  UCD-CCM-034,  Center  for  Computational 
Mathematics,  University  of  Colorado  at  Denver,  December  01  1994. 


69 


[12]  Achi  Brandt,  et.  al.  Algebraic  multigrid  (AMG)  for  sparse  matrix  equations. 
In  D.  J.  Evans,  editor.  Sparsity  and  Its  Applications,  Cambridge,  MA,  1984. 
Cambridge  University  Press. 

[13]  John.  W.  Ruge.  Algebraic  multigrid  (AMG)  for  geodetic  survey  problems.  In 
Steve  McCormick,  editor.  Preliminary  Proc.  Internal.  Multigrid  Conference,  Fort 
Collins,  CO,  1983.  Institute  for  Computational  Studies  at  Colorado  State  Univ. 

[14]  John.  W.  Ruge  and  K.  Stiiben.  Algebraic  multigrid.  In  Steve  McCormick,  editor. 
Multigrid  Methods,  pages  73-130,  Philadelphia,  Pennsylvania,  1987.  SIAM. 

[15]  Hans  Regler  and  Ulrich  Rude.  Layout  optimization  with  algebraic  multigrid 
methods  (AMG).  In  Proceedings  of  the  Sixth  Copper  Mountain  Conference  on 
Multigrid  Methods,  Copper  Mountain,  April  4-9,  1993,  Conference  Publication. 
NASA,  1993. 


70 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center . 2 

8725  John  J.  Kingman  Road.,  Ste  0944 

Ft.  Belvoir,  VA  22060-6218 

2.  Dudley  Knox  Library . 2 


Naval  Postgraduate  School 
411  Dyer  Rd. 

Monterey,  CA  93943-5101 

3.  Chairman,  Code  MA . 1 

Department  of  Mathematics 

Naval  Postgraduate  School 
1411  Cunningham  Road,  Rm  341 
Monterey,  Ca  93943-5216 

4.  Professor  Van  E.  Henson,  Code  MA/Hv . 2 

Department  of  Mathematics 

Naval  Postgraduate  School 
1411  Cunningham  Road,  Rm  341 
Monterey,  Ca  93943-5216 

5.  Lieutenant  Gerald  N.  Miranda,  Jr . 2 

Naval  Sumbarine  Base 

Submarine  Officer  Advanced  Course 
Code  N222,  SOAC  97060 
P.O.  Box  700 
Groton,  CT  06349-5700 


71 


